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I.  INTRODUCTION 


A.  CONCEPT 

The  purpose  of  this  thesis  is  to  characterize  typical  levels  of  vibrational 
noise  on  the  ocean  floor  to  ascertain  the  vibration’s  effect  on  possible  future 
bottom  mounted  sensors.  Lawrence  Livermore  National  Laboratory  (LLNL)  is 
interested  in  possibly  using  an  experimental  sensor  on  the  ocean  floor.  LLNL 
characterized  the  noise  response  for  various  terrestrial  versions  of  this  particular 
sensor,  and  they  expressed  interest  in  collaborating  in  a  study  of  the  vibrational 
noise  on  the  ocean  floor,  focusing  on  low  frequencies  during  quiescent  periods. 
The  data  used  for  this  thesis  was  obtained  from  publicly  available  recorded 
information  from  ocean  bottom  seismometers  (OBS).  It  is  beyond  the  scope  of 
this  thesis  to  characterize  all  low  frequency  noise  on  the  ocean  floor  since  noise 
can  differ  by  area  and  the  seismic  events  occurring  during  the  measurement.  The 
data  was  obtained,  then,  by  focusing  on  two  geographical  choke  points  which 
have  co-located  OBSs.  These  highly  trafficked  choke  points  were  considered  to 
be  a  good  representation  of  where  these  experimental  bottom  mounted  sensors 
might  be  located  should  they  be  built. 

Discussions  were  held  with  scientists  from  LLNL  to  find  out  sensor 
information  and,  specifically,  the  frequency  regimes,  amplitudes,  and  vibration 
axis  that  affect  the  current  sensors.  The  LLNL  personnel  discussed  what  an 
underwater  version  of  this  sensor  might  look  like,  but  they  were  more  interested 
in  the  statistics  of  the  vibrational  noise  of  the  ocean  floor  between  large 
amplitude  events  to  determine  what  to  account  for  in  future  designs. 

B.  RELATED  WORK 

A  thesis  by  Jonathon  P.  Scobo  entitled  “Ocean  Observing  Systems”  was 
used  for  this  research  to  start  obtaining  OBS  data  and  the  various  seismology 
software  to  process  it  [1].  Even  though  most  OBS  data  is  publicly  available,  OBS 
data  can  be  difficult  to  access  and  cumbersome  to  analyze.  Scobo’s  thesis 
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proved  essential  to  find  the  desired  OBSs,  which  seismology  computer  programs 
to  use  and  which  to  avoid,  and  how  to  request  specific  time  periods  of  data  [1]. 

Most  other  work  that  was  found  with  regards  to  seismic  vibrations  on  the 
ocean  floor  concentrated  on  larger  seismic  events  as  against  the  noise 
environment  between  these  events.  LLNL  expressed  less  interest  in  large  events 
and  more  interest  in  the  noise  environment  during  quiescent  periods  between 
large  amplitude  events. 

Somewhat  related  research  has  been  conducted  which  discusses  the 
shortcomings  of  different  types  of  OBSs  due  to  the  bottom  type  on  which  they  are 
mounted.  These  shortcomings  will  be  discussed,  because  the  data  used  for  this 
thesis  could  be  affected  by  these  shortcomings. 

C.  ORGANIZATION 

This  thesis  will  start  with  background  information  discussing  the 
Incorporated  Research  Institutions  for  Seismology  (IRIS)  website  and  the  Ocean 
Bottom  Seismograph  Instrument  Pool  (OBSIP).  The  background  chapter  also 
discusses  specific  OBS  sensors  used  for  this  thesis  and  OBS  fidelity  on  different 
bottom  types.  The  next  chapter  discusses  theory  (Chapter  III)  and  includes  an 
introduction  to  seismic  waves.  Chapter  IV,  Methodology,  discusses  the  process 
to  obtain,  process,  and  analyze  the  OBS  data.  Next,  the  results  chapter  (Chapter 
V)  contains  a  discussion  of  the  low  frequency  vibrational  noise  environment  in 
the  different  areas  of  interest.  Chapter  VI  discusses  the  overall  conclusions  as 
well  as  recommendations  on  future  work,  and  the  appendices  contain  the 
MATLAB  code  and  data  tables. 

D.  OVERVIEW 

Four  OBSs  in  two  different  geographic  locations  were  chosen  and 
subsequently  their  vibrational  noise  environments  were  investigated.  Figure  1 
from  the  IRIS  website  shows  a  global  view  of  publicly  available  OBSs.  This  global 
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view  allows  the  user  to  zoom  in  on  and  select  the  desired  OBS,  which  then 
opens  a  link  to  download  the  different  types  of  available  data. 


Figure  1.  The  Global  View  of  Publicly  Available  OBSs.  Source:  [2]. 

Two  geographic  choke  points  were  chosen:  the  Luzon  Strait  and  west  of 
the  Strait  of  Juan  de  Fuca.  IRIS  publicly  available  software,  specifically  a 
program  to  read  information  in  the  form  of  Standard  for  the  Exchange  of 
Earthquake  Data  (SEED),  called  “rdseed,”  a  program  used  to  evaluate  response 
files  called  “EVALRESP,”  and  Seismic  Analysis  Code  (SAC)  were  used  to 
transfer  downloaded  data  into  a  usable  format  and  then  calibrate  the  data.  The 
analysis  conducted  on  the  calibrated  data  was  accomplished  with  MATLAB. 
Figure  2  and  Figure  3  display  the  available  OBSs  at  the  chosen  choke  points. 
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Figure  2.  Available  OBSs  in  the  Luzon  Strait.  Source:  [2], 


Figure  3.  Available  OBSs  West  of  Strait  of  Juan  de  Fuca.  Source:  [2]. 

The  data  was  obtained,  then,  by  focusing  on  known  geographical  choke 

points  which  have  co-located  OBSs.  These  highly  trafficked  geographic  choke 

points  were  considered  to  be  a  good  representation  of  where  these  experimental 
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bottom  mounted  sensors  would  be  located  for  either  testing  or  operations  after 
they  are  built.  Initially,  more  geographic  choke  points  with  OBSs  were  inspected 
in  hopes  to  add  additional  locations  to  analyze.  One  of  the  OBS  locations  near  a 
choke  point  south  of  the  Bering  Strait  was  inspected,  but  the  OBSs  were  short 
period  OBSs  (discussed  later)  that  were  used  during  an  active  sound 
transmission  test.  At  this  location,  there  was  only  a  week’s  worth  of  noisy  data. 
According  to  IRIS  in  Figure  1  there  are  no  OBS  locations  for  the  majority  of  the 
operationally  significant  choke  points.  The  OBS  locations  are  understandably 
chosen  for  geological  significance  (e.g.,  hydrothermal  events  in  the  middle  of  the 
Atlantic  Ocean). 
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II.  BACKGROUND 


A.  IRIS  AND  OBSIP 

All  of  the  OBS  data  obtained  for  this  thesis  was  from  the  IRIS  and  OBSIP 
websites.  IRIS  is  a  non-profit  organization,  created  in  1984  by  the  National 
Science  Foundation  (NSF).  According  to  their  website,  IRIS  involves  over  100 
universities  in  the  United  States.  Their  programs  “contribute  to  scholarly 
research,  education,  earthquake  hazard  mitigation,  and  verification  of  the 
Comprehensive  Nuclear-Test-Ban  Treaty”  [3].  One  of  IRIS’s  mission  statements 
is  to  promote  seismic  data  availability,  and  they  allow  access  to  data  by  using  a 
simple  request  format,  which  will  be  addressed  later  in  this  thesis. 

IRIS  oversees  and  coordinates  multiple  seismic  networks,  which  includes 
the  network  used  for  this  thesis,  the  OBSIP.  OBSIP  focuses  on  marine 
seismology,  geology,  and  geodynamics.  The  OBSIP  website  [4]  also  contains 
links  describing  current  OBS  instrumentation,  links  to  access  OBS  data,  and 
information  for  researchers  to  request  OBSIP  instrument  usage.  The  OBSIP 
allows  their  OBS  equipment  to  be  available  to  NSF-sponsored  investigators  and 
to  other  government  research  or  educational  institutions  [4], 

B.  OBS 

Ocean  Bottom  Seismometers  can  be  broken  into  two  fundamental  types: 
long  period  and  short  period.  According  to  OBSIP’s  website  [4],  the  long  period 
OBSs  have  a  3-component  broadband  seismometer  (one  vertical  component  and 
two  horizontal)  and  a  low-frequency  differential  pressure  gauge 
(DPG)/hydrophone.  These  long  period  OBSs  are  primarily  designed  for  passive 
recording  and  are  normally  capable  of  a  greater  than  12-month  deployment.  The 
short  period  OBSs  consist  of  either  a  vertical  seismometer  or,  like  the  long  period 
OBS,  three  orthogonal  seismometers  and  a  hydrophone  [4],  The  OBSIP  website 
[4]  also  states  that  the  sampling  rates  for  the  short  period  OBSs  are  typically 
much  higher  (200-250Flz),  and  these  instruments  are  normally  deployed  from  1 
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to  6  months.  The  short  period  OBSs  are  typically  used  for  active-source 
experiments  [4], 

The  OBSs  used  in  this  thesis  are  long  period.  Because  of  how  OBSs  are 
placed  on  the  ocean  floor,  the  orientation  of  the  horizontal  channels  are  usually 
unknown.  Methods  to  orientate  the  horizontal  channels  can  involve  analyzing 
data  from  the  known  location  of  an  active  source  (e.g.,  air  gun)  or  using  P  or 
Rayleigh  waves  from  a  known  earthquake  [5].  The  orientation  of  the  horizontal 
channels  used  in  this  thesis  were  unknown. 

1.  LDEO  OBS  Sensor  MK2 

Figure  4  from  the  IRIS  website  displays  the  highlights  from  the  Lamont- 
Doherty  Earth  Observatory  (LDEO)  Standard  Ocean  Bottom  Seismometer.  The 
LDEO  OBSs  used  for  this  thesis  are  located  in  the  Luzon  Strait  and  were  used  in 
the  Taiwan  Integrated  Geodynamic  Research  (TAIGER)  project.  The  TAIGER 
project  is  a  joint  USA-Taiwan  program  that  studies  the  tectonic  development  of 
Taiwan  [6].  The  network  identifier  “YM”  designates  the  OBSs  in  the  TAIGER 
project.  The  two  OBS  stations  that  were  used  to  provide  data  were  YM  38  and 
YM  39,  displayed  in  Figure  5. 
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Ocean* Bottom  Seismology  Laboratory 

LDEO  Standard  Ocean-Bottom  Seismometer 


The  LDEO  standard  seismometer 
design  has  been  in  use  for  nearly  10 
years.  The  LDEO  OBS  lab  has  built 
and  operates  25  standard  OBSs  as 
part  of  the  NSF  OBS  Instrumentation 
Pool.  The  seismometer  sensor  is  an 
L4C  1  Hz  geophone,  with  a  low-noise 
amplifier,  giving  useful  response 
down  to  100  s,  and  a  differential 
pressure  gauge.  This  instrument  has 
been  used  in  both  year-long  passive- 
source  and  shorter-term  active- 
source  experiments.  The  design 
includes  dual  redundancy  with  two 
transponders  and  two  dropweights. 

Variants  and  add-ons 

Trawl-rcsistant  OBS 
Moored  DPG/hydrophonc 
Ocean-bottom  magnetometer 
Diffuse  flowmeter 
Trillium  Compact  sensor 
Absolute  pressure  gauge 
Hydrophone 


Specifications 

Max.  Depth 

5000  m 

Max.  Duration 

400  days  @  1 25  sps 

Channels 

4, 24-bit  recording 

Sensors 

L4C  3-component  geophones; 
differential  pressure  gauge 

Response 

1 00  s  -  60  Hz  (seismometer) 

0-20  Hz  (DPG) 

Leveling  system 

Active  360°,  motor-driven 

Weight 

750  lb  in  air 

Footprint 

3'X4' 

Flotation 

9X1 2”  glass  spheres 

Sampling 

40-100-125  sps 

Release 

Dual  dropweights 

Acoustics 

Two  ORE  1 2  kHz  transponders 

Power 

Lithium  battery  pack,  +/-  7.5  V 

Oscillator 

Seascan  10  MHz  clock 

Sensor  housing 

17’ glass  sphere 

Burnwires 

LDEO  design 

Recovery  aids 

Radio,  strobe,  flag 

Recording 

2  X  32  Gb  CompactFlash  cards 

Dropweights 

Two  steel  weights  (75  lb  in  air) 

Datalogger 

LDEO  ultra-low  power  OBS 
datalogger  (300  mW  @  1 25  sps) 

Deployment  History 


Figure  4.  Highlights  of  the  LDEO  OBS,  Which  Includes  a  Useful 
Instrument  Response  Down  to  a  Period  of  100  Seconds. 

Source:  [7], 
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Figure  5.  Location  of  YM38  and  YM39,  the  Two  OBS  Locations  in  the 

Luzon  Strait.  Source:  [2], 

The  LDEO  OBS  uses  an  L4C  1  Hz  geophone  as  the  seismometer  sensor 
[8],  which  is  displayed  in  a  drawing  in  Figure  6.  As  discussed  in  [9],  the  vertical 
L4C  geophone  produces  a  voltage  which  is  proportional  to  the  difference  in 
velocity  between  the  mass,  M,  and  the  ground  (x-xg).  The  mass,  M,  is 

suspended  with  a  spring  with  stiffness  k  and  a  dashpot,  b,  for  damping.  The 
voltage  is  produced  by  the  motion  sensing  coil  with  inductance,  Lc,  and  coil 
resistance,  Rc.  The  L4C  1  Hz  geophone  has  a  harmonic  frequency  of  1  Hz,  and 
the  resistor,  Rs,  is  an  external  resistor  to  shunt  the  output  and  critically  damp  the 
system  at  that  frequency,  as  explained  in  the  paper  by  Bowden  [9], 
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Figure  6.  Drawing  of  Vertical  L4C  1  Hz  Geophone,  Which  Is  Used  as  the 
Seismometer  Sensor  of  the  LDEO  OBS.  Source:  [9]. 

2.  SIO  ABALONES  OBS  System 

The  Scripps  Institution  of  Oceanography  (SIO)  Autonomous  Broad 
Application  Low  Obstruction  Noise  Exempt  System  (ABALONES)  OBS  system 
can  have  a  variety  of  configurations.  The  particular  OBSs  used  for  this  thesis, 
located  to  the  west  of  the  Strait  of  Juan  de  Fuca  as  displayed  in  Figure  7,  are 
considered  an  intermediate-period  sensor  because  of  their  large  frequency 
range.  In  this  configuration  they  concentrate  on  long-period  vibrations  with  a 
Nyquist  frequency  of  25  Hz.  Table  1  displays  the  specification  for  the  SIO 
ABALONES  system,  and  Figure  8  displays  a  3D-CAD  drawing  of  the  ABALONES 
OBS  design.  The  yellow  shell  of  the  ABALONES  design  is  a  trawl  and  current 
resistant  enclosure,  allowing  these  OBSs  to  be  set  in  shallow,  continental  shelf 
areas  to  as  deep  as  6000m  [10].  It  is  equipped  with  a  Nanometrics  Trillium 
Compact  Seismometer,  which  has  the  advantage  of  a  nearly  flat  response  curve 
from  1 20  seconds  to  1 00  Hz  [1 0],  [1 1  ]. 
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Figure  7.  7D  J73A  and  7D  M02A,  the  Two  OBS  Locations  West  of  the 
Strait  of  Juan  de  Fuca  Used  for  this  Thesis.  Source:  [2], 


Table  1 .  SIO  ABALONES  OBS  Specifications.  Source:  [11]. 


Seismometer/ 

Accelerometer 

Intermediate-Period:  Nanometrics  Trillium-Compact  seismometer,  which 
has  a  flat  velocity  response  from  120s  to  100  Hz.  The  seismometer  is 
housed  within  an  8”  Inside  Diameter  (ID)  Titanium  cylinder  and  a  custom- 
designed  active  leveling  system  which  can  be  leveled  from  ±180°  to 
vertical  with  ~0.3°  accuracy. 

Pressure  Sensor 

Short-period:  Customized  HiTech  hydrophone  with  internal  preamp. 
Bandwidth  (-3  dB)  is  50  mHz  -  15  kHz. 

Long-period:  Differential  Pressure  Gauge  (DPG)  with  response  from  10 
mHz  -  10  Hz. 

Digitizer 

24-bit  A/D  with  solid-state  recording  (CompactFlash).  The  dynamic  range 
is  ~  126  dB  (3-bit  self-noise).  The  ABALONES  design  includes  the  latest 
seismic  A/D  chip  with  programmable  sample  rates  from  1  Hz  to  4  kHz. 

Clock 

Low-power,  digitally-temperature-compensated  (DTCXO),  precision  time 
base.  The  drift  rate  is  1:3-5  x  10  ®  (<5  ms/day  before  correction  and  <1.5 
s/yr) 

Recording 

Capacity 

Single  CompactFlash  card  slot  ->  64  Gb  CF  cards  currently  available. 
Capable  of  recording  1-year  deployment  on  4-channels  @  100  Hz. 

Data  Offload 

USB2  through  the  end  cap  without  opening  the  pressure  case.  Offload 
rates  are  CF  card  dependent,  currently  up  to  400  Mb/s. 

Battery  Pack 

Batteries  are  mounted  in  the  main  data  logger  pressure  case  and  optionally 
in  an  additional  battery  case.  Alkaline  batteries  power  the  acoustic  release. 

Recording 

duration 

IPA:  Lithium  battery  pack  can  provide  power  for  12  months,  Alkaline 
packs  can  provide  power  for  4+  months. 

SPA:  Lithium  battery  pack  can  provide  power  for  18+  months. 

Weight 

850/550  lbs  with/without  anchor  (in  air) 

Pressure  Case 

5”  diameter  A1  cylinder  -  6  km  depth  rating 

Release 

Double  burn  wire  operated  acoustically 

Dimensions 

29"  high  x  42"  x  42"  (w/  bail) 

Power  (Total) 

Intermediate:  <300  mW  (4  chan),  plus  ~160  mW  w/  Trillium-Compact 
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Figure  8.  3D-CAD  Drawing  of  ABALONES  Ocean  Bottom 
Seismometer  Design.  Source:  [11]. 


C.  OBS  FIDELITY 

As  mentioned  by  Sutton  et  al.,  “The  usual  reason  for  placing  a 
seismometer  anywhere  is  to  measure  free-field  particle  motion  in  response  to 
some  seismic  source”  [12].  The  complication  presented  by  soft  sediment  of  the 
ocean  floor  is  that  the  coupling  between  the  seismometer  and  the  sediment  may 
be  unknown  [12],  Vertical  vibration,  channel  Z  in  this  thesis,  is  of  more  concern 
than  any  horizontal  motion,  channels  1  and  2.  Duennebier  and  Sutton  maintain 
that  vertical  motion  across  the  water-sediment  boundary  is  reasonably 
represented  by  an  OBS  due  to  particle  motion  continuity  across  the  interface 
[12],  [13].  The  fidelity  of  horizontal  motion  captured  by  an  OBS  is  also  less 
accurate  because  of  the  shear  discontinuity  between  the  sediment  and  the  water 
[14].  Problems  with  OBS  horizontal  fidelity  are  compounded  by  bottom  ocean 
currents  in  addition  to  the  soft  sediment. 
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In  [14],  Duennebier  and  Sutton  point  out  that  horizontal  motion  is 
discontinuous  across  the  boundary  of  the  ocean  floor.  This  differential  motion 
generates  torque  on  the  instrument  which  in  turn  causes  both  horizontal  and 
vertical  motion  to  be  captured  by  the  seismometer  channels.  In  addition,  if  the 
sensor  is  not  completely  level,  components  of  vertical  and  horizontal  motions  will 
couple  into  the  wrong  channels.  Consequently,  cross-correlation  between  the 
horizontal  and  vertical  channels  is  common  on  bottom  mounted  OBSs  [14].  The 
horizontal  channels  are  considered  to  be  more  susceptible  to  noise  from  this 
source.  Duennebier  and  Sutton  collected  data  in  a  study  proving  increased  OBS 
fidelity  when  the  OBS  is  buried  beneath  the  sediment  on  the  ocean  floor  [13]. 

Any  sensor  placed  on  the  ocean  floor  will  be  subject  to  the  same 
vibrations  as  the  OBS.  Therefore,  these  issues  relating  to  the  fidelity  of  the 
seismic  measurements  may  be  of  concern  to  LLNL  for  their  potential  to  affect  the 
performance  of  the  proposed  sensor  as  well  as  for  possible  limitations  in  the 
accuracy  of  the  seismic  measurements  provided  in  this  thesis.  If  a  bottom- 
mounted  sensor  is  used,  [14]  explains  that  the  effects  of  unwanted  vertical 
motion  can  be  minimized  by  increasing  the  footprint  of  a  sensor  package  relative 
to  its  height.  The  tradeoff  is  a  decreased  sensitivity  to  certain  types  of  waves 
when  the  sensor  base  approaches  too  high  a  fraction  of  the  seismic  wavelength, 
particularly  in  higher  frequency  shear  waves. 
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III.  THEORY 


Seismic  waves  can  be  divided  into  two  categories:  body  waves  and 
surface  waves.  The  body  waves  travel  through  the  interior  of  the  earth  and  are 
the  first  waves  to  arrive  after  a  seismic  event  [15].  The  surface  waves  cause  the 
most  damage  associated  with  an  earthquake,  and  they  travel  along  the  earth’s 
crust  [15].  The  mathematical  theory  behind  both  types  of  waves  is  discussed  in 
this  chapter.  Figure  9  is  the  coordinate  system  used  in  the  equations.  All  of  the 
equations  in  this  chapter  are  from  Bullen  and  Bolt  [16]. 


*3 


*1 


Figure  9.  Coordinate  System  Used  in  Equations  1  through  20. 


A.  BODY  WAVES 
1 .  P  Wave 

The  P  wave  is  the  primary  or  “push  wave.”  Bullen  and  Bolt  [16]  describe 
the  P  wave  in  terms  of  the  dilatation,  9 .  The  dilatation,  or  irrotational  disturbance, 
is  the  fractional  change  in  volume  of  a  point  in  the  medium  as  the  wave  passes.  It 
is  given  by  the  divergence  of  the  vector  ui ,  which  gives  the  displacement  of  a 

point  with  respect  to  its  equilibrium  position  in  direction  xi  as  shown  in  Equation 

1 .  Equation  2  is  the  scalar  wave  equation  for  the  dilatational  disturbance.  The  P 
wave  is  transmitted  through  a  substance  at  speed  a  ,  given  in  Equation  3. 


_  du,  du ,  du, 

9  =  — L  +  — -  + — 1 
dxx  dx2  dx3 


(1) 
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(2) 


d2e 

dt2 


a2y26 


a  = 


(  A  +  2//V2 

P  j 


(3) 


A  is  Lame’s  first  parameter,  which  is  an  elastic  modulus  defined  by 
Equation  4.  p  is  Lame’s  second  parameter,  which  is  also  known  as  the  shear 
modulus,  and  p  is  the  volume  density  of  the  material.  The  Lame  parameters 
together  characterize  the  linear  elasticity  for  homogeneous  isotropic  media  [17]. 
The  relationship  between  Lame’s  first  parameter  and  the  shear  modulus,  p ,  is 
given  by  [16] 

A  =  K~p  =  p(a2 -2/32),  (4) 

where  K  is  the  bulk  modulus  of  the  material. 


Figure  10  illustrates  how  a  P  wave  travels.  It  is  a  longitudinal  wave 
consisting  of  changes  in  the  density  of  the  material. 


T  =  0 


T  =  1 


T  =  2 


T  =  3 


Figure  10.  Illustration  of  P  Wave.  Source:  [15]. 
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2. 


S  Wave 


S  waves  are  a  type  of  shear  wave  and  are  described  by  the  vector  wave 
equation 


dr 


curl  (u ,)  =  /32V2curl(uj) 


(5) 


where  J3  is  the  speed  of  the  rotational  disturbance  [16]  and  given  by 


fi  = 


'Ey 

Ip; 


(6) 


Since  the  speed  depends  on  the  shear  modulus,  S  waves  will  not 
propagate  through  fluids  which  cannot  support  shear  stress.  The  S  wave  in  the 
far  field  can  be  plane  polarized  [16].  If  an  S  wave  is  traveling  horizontally  and  the 
particle  motion  is  also  horizontal,  it  is  denoted  SH.  If  the  particle  motion  is  vertical 
only,  the  S  wave  is  denoted  SV  [16].  Figure  11  illustrates  an  SV  wave  assuming 
that  z  is  vertical  and  x  and  y  are  horizontal. 
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T  =  0 

T  =  1 

T  =  2 

T  =  3 


Figure  1 1 .  Illustration  of  an  S  Wave.  Source:  [1 5]. 

B.  SURFACE  WAVES 

In  Figure  12,  which  is  adapted  from  an  illustration  in  Bullen  and  Bolt  [16], 
M  is  an  elastic,  homogeneous  half-space  below  the  plane  denoted  by  the  x1  and 

x,  axes.  Mis  bordered  above  the  xl/x2  plane  by  another  homogeneous  elastic 
half-space,  M’. 


Figure  12.  Two  Isotropic  Homogeneous  Elastic  Media  M’  and  M. 

Adapted  from  [16]. 
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Waves  can  travel  along  the  boundary  between  these  two  half  spaces.  To 
be  trapped  along  the  boundary,  the  amplitude  of  these  waves  must  decrease 
exponentially  as  the  distance  from  the  boundary  increases.  Consider  such  a 
plane  wave  traveling  in  the  x1  direction.  Since  it  is  a  plane  wave,  the  particle 

displacement  along  any  line  parallel  to  x2  is  equal.  As  a  result,  any  partial 
derivative  taken  with  respect  to  x,  will  be  zero  [16].  Letting  <j)  and  y/  be  functions 


which  can  be  used  to  express  the  displacement  components  ux  and  u3  as  given 
by 


d(f>  By/  d(/>  dy/ 

/'/j  I  ^  IA> 3  j 

dxx  dx3  dx3  dxx 


(7) 


it  can  be  shown  that 


V2(j)  =  0 
du ,  du3 

v  V  =  — 1 — 1 

dx3  dx] 


(8) 

0) 


This  suggests  that  ^5  is  linked  with  P  waves  since  they  are  described  by 
the  dilatation,  y/  is  linked  with  SV  waves  since  it  is  associated  with  a  differential 
displacement  between  the  direction  of  propagation  and  the  vertical.  Any 
displacement,  u2,  would  be  linked  with  SH  waves. 


The  wave  equations  associated  with  <j) ,  y/  ,  and  u2  are  given  by 


d2(j) 
8t 2 


=  a2V2<f> 


82y / 

~gF 

d2u2 

~w 


=  J32V V 

=  p2v\ 


(10) 

(11) 

(12) 


where  the  speeds  in  M’  are  replaced  with  a'  and  /?'. 
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To  find  solutions  for  y/  ,  and  u2  [16]  sets  r  and  5  to  the  following, 

where  c  is  the  speed  of  a  set  of  simple  harmonic  waves  with  wavelength  2tt I  k  . 
Again,  for  M\  the  variables  are  r'  and  s'.  The  variables  r,  r' ,  s,  and  5' 
designated  in  Equations  13  and  14  are  imaginary. 


r  = 


f  2 


c 


^\2 

1  ,5  = 

) 


f  2 


C 


1  1 


(  2  X 

c 

2 

f  2  X 

c 

- T-l 

,S  ’  = 

— T-l 

U’2  J 

[p'2  j 

Solutions  in  M  and  M’  are  sought  in  the  form: 

(j)  _  AgikC-rxs+^-ct) 

y/  =  Be'k(~sx?,+x'~c  r) 

U  =  (Jgik('-SX3+XL~Ct) 

Q  —  Q'  e‘k(-r'x3+x,-ct) 

yj  —  J?'  giki.-s'ts+Xi-ct) 

u  —  p'  eik{-s'xl+xl-ct) 


(13) 

(14) 

(15) 

(16) 

(17) 

(18) 

(19) 

(20) 


Since  both  r  and  5  are  imaginary,  the  positive  or  negative  square  root  is 
taken  as  needed  to  ensure  that  all  solutions  consist  of  waves  that  exponentially 
decrease  as  the  distance  away  from  the  boundary  is  increased.  Expressions  for 
the  displacements  ut  and  u3  as  a  function  of  time  and  position  can  be  found  from 
these  solutions  by  applying  the  relationships  which  implicitly  define  ^  and  y/  in 
Equation  7. 

The  boundary  conditions  dictate  which  solutions  propagate.  For  example, 
if  M  is  not  homogeneous  elastic  half  space,  but  is  instead  a  vacuum  (or  air)  the 
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resulting  wave  is  a  Rayleigh  wave,  where  there  are  no  SH  waves.  The  particle 
motion  is  retrograde,  and  displayed  in  Figure  13. 


T  =  0 


T  =  1 


T  =  2 


T  =  3 


Figure  13.  Illustration  of  Rayleigh  Wave.  Source:  [15]. 

If  M  is  not  a  homogeneous  elastic  half  space,  but  is  instead  an 
incompressible  fluid  (water)  the  resulting  wave  is  a  Scholte  wave  [18].  The 
particle  motion  again  is  retrograde,  and  the  amplitude  decays  exponentially  with 
distance  from  the  boundary.  An  illustration  of  a  Scholte  wave  is  displayed  in 
Figure  14. 


Figure  14.  Illustration  of  Scholte  Wave.  Source:  [19]. 
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IV.  METHODOLOGY 


A.  OVERVIEW 

This  chapter  explains  the  methods  chosen  to  calibrate  the  data  using  IRIS 
software,  why  these  methods  were  chosen,  and  the  data  analysis  conducted  with 
MATLAB.  Multiple  routes  to  calibration  were  available,  and  picking  the  most 
effective  way  forward  proved  challenging.  Many  different  software  calibrations 
methods  are  available  from  IRIS,  and  a  few  different  methods  to  get  calibrated 
data  on  MATLAB  (manually)  were  attempted  before  the  steps  were  decided 
upon.  The  steps  laid  out  in  this  chapter  should  allow  follow-on  study  with 
somewhat  less  of  a  learning  curve  for  non-seismologists. 

B.  CHOOSING  THE  OBS 

The  global  view  on  the  IRIS  website,  displayed  in  Figure  1,  of  the  available 
OBSs  proved  a  convenient  place  to  start  the  OBS  selection  process.  This  global 
view  allows  the  user  to  focus  on  the  desired  area  and  select  individual  sensors. 
Once  selected,  a  window,  which  is  shown  in  Figure  15,  is  opened  that  displays  to 
the  user  the  dates  of  available  data  and  a  link  to  more  information.  It  becomes 
readily  apparent  if  the  selected  OBS  is  a  short  period  vice  long  period.  In  most 
cases,  the  short  period  OBS  dates  available  will  be  for  days  and  weeks  versus 
multiple  months. 
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Figure  1 5.  Convenient  Use  of  Global  View  for  OBS  Selection. 

Source:  [1]. 

Once  the  user  clicks  on  the  “more  information”  link,  the  specific  OBS  page 
opens  and  allows  the  user  to  view  the  network  and  station  information  with  the 
channels  available.  Figure  16  shows  this  page.  The  user  can  select  any  of  the 
channels  and  view  the  response  curve,  the  depth,  sensor  location  in  latitude  and 
longitude,  and  the  response  (RESP)  and  poles  and  zeroes  (PZ)  files. 
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JRIS^  IRIS  DMC  MetaData  Aggregator 

Station  summary  (1  time  span) 


Legend:  QQQQ 


Network 

YM  ::  TAIwan  Integrated  GEodvnamic  Research  (TAIGERl ::  YM  Network  Map  ::  DOI 

Station 

39  ::  LDEO  OBS  39  ::  TAIGER  ::  39  Station  Map  ::  RESP  ::  SAC  PZs  ::  XML 

Latitude 

21.599700 

Longitude 

121.364100 

Elevation 

-2508 

Start 

2008/05/11  (132)00:00:00 

End 

2009/05/30  (150)  23:59:59 

Epoch 

2008/05/11  (132)  21:25:51  -  2009/05/30  (150)  21:20:11 

Instrument 

LDEO  OBS  Sensor  Mk2/LDEO_new  dummy  at  40/125 

Channels  (Hz) 

Location  BHJ.  (40)Q,  BH2  (40)Q,  BHZ  (40)0 

Instrument 

LDEO  DPG/LDEO_new  dummy  at  40/125 

Channels  (Hz) 

Location  BDH  (40)Q 

MetaData  Load 

2015/03/30(089)  14:25:20 

Virtual  network  affiliations: 

Name  Description  Primary  DC  Secondary  DC 

OBSIP  U.S.  National  Ocean  Bottom  Seismograph  Instruments  OBSIP  IRIS  DMC 

■OBSIP  LDEO  OBSIP  Experiments  from  LDEO  OBSIP  IRIS  DMC 

OBSIP  LDEO  OBSIP  Experiments  from  LDEO  OBSIP  IRIS  DMC 

.UNRESTRICTED  All  unrestricted  stations,  generated  via  cron  IRIS  DMC  IRIS  DMC 

No  data  available  in  real-time  systems 

Archive  data  availability  -  Make  a  batch  request  for  data  fbreq  fasti  -  (data  access  overview) 

Earliest  Latest 

O  2008/05/11  (132)22:25:51  2009/05/30(150)21:20:11 

Query  archive  for  data  day  availability: 

2008  2009 


Figure  16.  Detailed  Information  on  OBS  Selected  from 

Global  Map.  Source:  [20]. 

C.  REQUESTING  DATA  FROM  IRIS 

There  are  multiple  ways  to  request  data  from  IRIS.  The  method  used  in 
this  thesis  made  use  of  the  link  in  Figure  16,  “Make  a  batch  request  for  data 
(breq_fast).”  The  breq_fast  request  form,  displayed  in  Figure  17,  is  self- 
explanatory  with  required  items  that  must  be  entered  starred.  The  “label”  can  be 
anything  the  user  wants,  which  is  the  name  of  the  file  sent  by  the  Data 
Management  Center  (DMC)  to  the  entered  email.  The  DMC  is  the  part  of  IRIS 
that  stores  and  makes  available  the  seismic  data  [21].  After  the  “Start  Query” 
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button  is  pressed,  Figure  18  depicts  the  final  webpage  opened.  For  this  thesis, 
“Full  SEED”  was  the  format  used.  After  the  form  is  submitted,  IRIS  sends  an 
email  with  a  link  to  download  the  requested  SEED  file  directly.  The  link  is 
available  for  24  hours.  “Full  SEED”  is  used  because  of  the  various  calibration 
files  that  Full  SEED  makes  available  with  IRIS  software. 


summaries 

by  station 

by  network 

by  timeseries 

virtual  nets 

breq_fast 

channels 

stations 

responses 

temp  networks 

assembled 

events  comments 

BREQ_FAST  Request  Form 


virtual  network 

network  ym 


latitude  and  longitude 
NORTH 


station  39 
location 
channel 

data  start  time*  2008  Jun  Q  13  3  000000 

data  end  time*  2008  Jun  @  13  3  235959 


sample  rate  >=  and  <= 
flags  like 
comments  like 
sensor  type  like 
site  like 

data  quality  Best 


channel  parameters 


WEST 

u 

EAST 

SOUTH 


fcjearH 

elevation 

depth 

>= 

>= 

<= 

<= 

azimuth 

dip 

>= 

>= 

<= 

<= 

name*:  Jeremy  Hankins 
institution*:  Naval  Postgraduate  School 
Street  address*:  I  University  Cir.,  Monterey, 
email*:  I 
phone*: 
label*: 

media:  Electronic  (FTP)  @ 

alternate  media:  dat  4mm 

alternate  media:  dit  jjjjJ 

DMC  archived  waveforms  or 
query  over  (  metadata/dataiess/RESP«* 

Start  Query 

•These  fields  arc  required. 

♦*A  query  on  metadata  may  result  in  a  request  file  for  unavailable  data,  but  returns  much  faster. 


Figure  17.  The  BREQ_FAST  Request  Form  Used  to  Get  SEED  Data 

from  IRIS.  Source:  [22], 
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Notice 


If  this  query  takes  mote  than  5  minutes,  you  may  wish  to  hit  the  STOP  button  on  this  browser  window,  further  restrict  your  query,  and  resubmit  the 
query. 

Querying  network  like  'YM'  for  year=2008 
Found  waveform  data  for  1  stations  and  4  channels. 


Make  Map  of  Stations 


Submit  request  below  to  receive 


full  SEED 
datalcss  SEED 
mini  SEED 
sync  file 
RESP  file 

email  request  to  myself 


.NAME  Jeremy  Hankins 

•INST  Naval  Postgraduate  School 

•MAIL  1  University  Cir.,  Monterey,  CA  93943 

.EMAIL  jhankinsPnps.edu 

.PHONE  9045044708 

•FAX  YOUR_FAX 

•MEDIA:  FTP 

.ALTERNATE  MEDIA:  DAT 

.ALTERNATE  MEDIA:  DLT3 

.LABEL  YM_39_13 JUN08 

.FROM  SO 

.QUALITY  B 

.END 


39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BDH 

(View. 

in 

MDA). 

39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BH1 

(View 

in 

MDA). 

39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BH2 

fVicw 

in 

MDA1 

39 

YM 

2008 

06 

13 

00 

00 

00.0000 

2008 

06 

13 

23 

59 

59.9999 

1 

BHZ 

fView 

in 

MDA1 

Figure  18.  Final  Submission  Page  to  Get  SEED  Data  from  IRIS. 

Source:  [23]. 


D.  DATA  ANALYSIS 

It  is  at  this  point  where  a  Mac  computer  (or  another  Unix-based  machine) 
is  required  to  process  the  data.  As  pointed  out  in  Scobo’s  thesis  [1],  Unix  is 
required  to  run  IRIS’s  available  software.  A  Mac  computer  proves  useful  because 
of  the  UNIX  availability  in  the  “Terminal”  application.  With  an  Internet  search  of 
Unix  commands  combined  with  the  manuals  that  are  available  on  IRIS’s  website, 
a  novice  Unix  user  can  at  least  use  the  installed  programs  from  IRIS  with 
practice.  Installing  and  compiling  Unix  code,  though,  is  not  a  trivial  endeavor 
without  previous  Unix  experience,  and  it  is  recommended  to  get  assistance 
installing  and  compiling  from  someone  with  skills  in  Unix  before  proceeding  with 
the  data  analysis  portion.  The  IRIS  website  conveniently  has  instruction  manuals 
for  all  of  its  software  with  examples  that  greatly  assist  novice  operators. 
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1. 


Rdseed  v5.3.1 


Rdseed  reads  the  SEED  files  and  outputs  the  raw  data  in  SAC  format, 
which  is  a  binary  format  that  can  be  read  by  the  SAC  software  or  MATLAB.  In 
addition,  RESP  files  and  PZ  files  for  all  available  channels  can  be  output  by 
rdseed  for  use  in  the  calibration  process.  Rdseed  proved  effective  and  somewhat 
user-friendly  after  installed.  Figure  19  displays  a  terminal  window  with  typical 
rdseed  commands.  It  is  recommended  to  have  an  interactive  manual  open  during 
the  first  few  times  using  this  software.  The  user  can  choose  the  dates  and  times 
to  limit  the  SAC  file’s  size.  A  user  could  also  download  a  multiple-day  SEED  file 
and  then  use  rdseed  to  chop  the  SEED  file  into  the  desired  times  for  SAC  file 
output.  The  rdseed  depicted  in  Figure  19  is  in  the  prompt  mode,  which  walks  the 
user  through  the  inputs  to  get  a  SAC  file,  a  RESP  file,  and  a  PZ  file  (if  desired). 


Iphs-imac: rdseedv5.3.1  phimac7$  ./rdseed 
«  IRIS  SEED  Reader,  Release  5.3.1  » 

Input  File  (/dev/nrst0)  or  'Quit*  to  Exit:  /users/phimac7/desktop/YM_39/YM_39_DEC_08/YM_39_DE 
C_13_14. 104940. seed 
Output  File  (stdout) 

Volume  #  [ < 1 ) — N ] 

Options  lacCsSpRtde] 

Summary  file  (None) 

Station  List  (ALL) 

Channel  List  (ALL) 

Network  List  (ALL) 

Loc  Ids  (ALL  (" — "  for  spaces])  : 

Output  Format  [(1=SAC),  2=AH,  3=CSS,  4=mini  seed,  5=seed,  6=SAC  ASCII,  7=SEGY,  8=Simple  ASCII( 

SLIST),  9=Simple  ASCII(TSPAIR) ]  :  1 

Output  file  names  include  endtime?  [Y/(N)]Y 

Output  poles  &  zeroes  ?  [Y/(N)]N 

Check  Reversal  [(0=No),  l=Dip. Azimuth,  2=Gain,  3=Both] :  0 

Select  Data  Type  [(E=Everything) ,  D=Data  of  Undetermined  State,  M=Merged  data,  R=Raw  waveform 
Data,  Q=QC'd  data]  :E 

Start  Time(s)  YYYY, DDD, HH:MM: SS. FFFF  :  2008,348,00:00:00 
End  Time(s)  YYYY, DDD, HH:MM: SS. FFFF  :  2008,348,23:59:59.9999 
Sample  Buffer  Length  [20000000]: 

Extract  Responses  [Y/(N)J  :  Y 

Output  data  format  will  be  sac. binary. 

INFO:  sac  variable  EVDP  will  be  in  KILOMETERS 

Writing  YM.39..BDH,  3456000  samples  (binary),  starting  2008,348  00:00:00.0080  UT 

Writing  YM.39..BH1,  3456000  samples  (binary),  starting  2008,348  00:00:00.0080  UT 

Writing  YM.39..BH2,  3456000  samples  (binary),  starting  2008,348  00:00:00.0080  UT 

Warning...  Azimuth/Dip  Reversal  found  39.BHZ,  Data  inversion  was  not  selected 

Writing  YM.39..BHZ,  3456000  samples  (binary),  starting  2008,348  00:00:00.0080  UT 

Writing  RESPONSE  file:  RESP. YM. 39. . BOH 
Writing  RESPONSE  file:  RESP.YM.39. .BH1 

.  BH2 
.  8HZ 


Writing  RESPONSE  file:  RESP.YM.39.. 
Writing  RESPONSE  file:  RESP.YM.39. 


Input  File  (/dev/nrst0)  or  'Quit'  to  Exit:  [] 


Figure  19.  Example  Rdseed  Usage  in  Prompt  Mode. 
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2. 


EVALRESP 


For  this  thesis,  EVALRESP  was  used  to  produce  frequency  amplitude 
phase  (FAP)  files.  The  FAP  file  is  an  ASCII  file  with  three  columns  (frequency, 
amplitude,  and  phase)  which  can  be  used  to  apply  the  sensor  response  and 
obtain  calibrated  data.  EVALRESP  is  built  into  the  SAC  program’s  “transfer” 
command,  but  if  the  user  wants  to  obtain  a  FAP  file,  the  EVALRESP  software 
must  be  run  separately.  According  to  the  SAC  Manual  [24],  the  advantage  of 
using  a  FAP  with  the  “transfer”  function  “is  that  one  can  include  additional  stages 
of  the  instrument  response  and/or  control  more  explicitly  the  frequency  range 
over  which  the  correction  is  applied.”  A  frequency  range  is  chosen  when  building 
a  FAP  file  within  EVALRESP. 

Figure  20  depicts  a  sample  EVALRESP  command  to  generate  a  FAP  file 
from  the  YM  39  OBS  used  in  this  thesis.  Because  channels  1,  2,  and  Z  were 
used  in  data  analysis,  a  FAP  file  for  all  three  channels  had  to  be  generated  in  the 
method  shown  in  Figure  20  for  channel  Z.  The  EVALRESP  manual  [25]  provides 
detailed  information  on  the  commands  in  Figure  20.  In  this  figure,  “39”  is  the 
station,  “BHZ”  is  the  channel,  2008  is  the  year,  348  is  the  Julian  date,  “0  20’  lists 
the  frequency  range  sought  after,  “2049”  is  the  requested  number  of  samples  in 
the  given  frequency  range,  and  “-f”  precedes  the  file  location.  The  RESP  file 
called,  in  this  case  “RESP.YM.39..BHZ,”  is  the  response  file  generated  by  the 
rdseed  program,  “-s”  is  the  type  of  spacing  requested  (“tin”  for  linear  spacing  is 
chosen  here),  and  “-r”  is  the  response  type  which  is  set  to  “FAP”  [25]. 

phs-imac:src  phimac7$  ./evalresp  39  BHZ  2008  348  0  20  2049  -f  /Users/phimac7/desktop/YM_39/YM_39 
_DEC_08/RESP.YM.39..BHZ  -s  tin  -r  'tap' 
phs-imac:src  phimac7$  0 

Figure  20.  Sample  EVALRESP  Input  for  FAP  File. 

In  addition  to  a  FAP  file,  EVALRESP  can  also  output  a  complex  spectra 
(CS)  file  and  an  amplitude  phase  (AP)  file.  The  CS  file  is  a  three  column  file 
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containing  the  frequency,  real  part  of  the  amplitude,  and  imaginary  part  of  the 
amplitude,  whereas  the  AP  file  is  similar  to  the  FAP  file  except  that  frequency  is 
not  included. 

3.  SAC 

SAC  software,  with  the  built-in  transfer  function,  produced  the  calibrated 
data  for  this  thesis.  Figure  21  demonstrates  the  commands  used  for  this  thesis. 
Again,  it  is  recommended  to  open  an  interactive  software  manual,  available  on 
the  IRIS  website,  when  using  SAC.  After  reading-in  the  SAC  data  file  with  the 
read  (r)  command,  any  linear  trends  and  offsets  were  removed  with  a  remove 
trend  command  (rtr).  The  digital  Fourier  transform  (DFT)  that  is  used  with  a 
transfer  function  can  pad  the  beginning  and  end  of  a  time  series  with  zeros  so 
that  SAC  can  achieve  a  power  of  2  number  of  points.  The  “taper”  command  was 
applied  to  the  data  sequence  after  removing  the  trend  to  reduce  sidelobes  in  the 
spectral  data.  The  default  setting  for  the  taper  command  is  a  Hanning  taper  with 
a  width  of  .05  times  the  entire  length  of  data  used  [24],  Therefore,  the  amplitude 
of  10%  of  the  data  sequence  is  reduced  by  the  taper  and  was  discarded  from 
analysis. 

The  SAC  program  has  multiple  built-in  seismic  analysis  tools. 
Unfortunately,  the  SAC  manual  does  not  go  into  much  detail  discussing  the 
theory  behind  how  these  analysis  tools  operate.  The  “transfer”  command  built 
into  the  SAC  program  “performs  deconvolution  to  remove  an  instrument 
response  and  convolution  to  apply  another  instrument  response.”  [24]  The 
transfer  command  can  be  achieved  using  either  a  RESP  file,  a  FAP  file,  or  a  PZ 
file.  During  a  transfer  with  FAP  file  option,  if  a  frequency  is  encountered  below 
the  lowest  frequency  in  the  FAP  file,  the  algorithm  uses  the  lowest  frequency 
given  for  the  correction.  Similarly,  if  the  frequency  encountered  is  above  the  max 
FAP  file  frequency,  the  algorithm  uses  the  highest  frequency  given.  The  rdseed 
program  produces  the  PZ  and  RESP  files.  The  EVALRESP  program,  as 
discussed  earlier,  produces  the  FAP  file.  The  four  numbers  following  “freq”  in 
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Figure  21  are  the  frequency  limits,  fl  through  f4.  According  to  the  SAC  Manual 
[24],  the  frequency  limits  should  be  such  that  fl  <  f2  <  f3  <  f4.  A  high-pass  taper 
is  applied  between  fl  and  f2,  no  taper  (unity)  is  applied  between  f2  and  f3,  and  a 
low-pass  taper  is  applied  between  f3  and  f4.  To  avoid  ringing,  [24]  recommends 
that  f2=2*f1  and  f4>=2*f3.  When  using  SAC,  the  file  that  is  read-in  will  remain 
there  unless  the  program  is  commanded  to  do  otherwise.  The  INICM  command 
deletes  any  data  in  memory  and  was  used  prior  to  any  new  file  to  be  read  in  [24], 


SEISMIC  ANALYSIS  CODE  [11/11/2013  (Version  101.6a)) 

Copyright  1995  Regents  of  the  University  of  California 

SAC>  r  /users/phimac7/desktop/7D_J65A/NOV2011/NOV15/2011. 319.00.00. 00. 0195_2011. 319. 2 
3. 59. 59. 0195. 7D. J65A. . BHZ.M. SAC 
SAC>  rtr 
SAC>  taper 

SAC>  trans  from  fap  s  /users/phimac7/desktop/FAP. 70. J65A. .BHZ  to  vel  freq  .01  .02  10 
20 

Extrapolating  below  lowest  FAPFILE  frequency 
Extrapolating  above  highest  FAPFILE  frequency 
Station  (J65A  ),  Channel  (BHZ  ) 

SAC>  write  NOV15_2011_7D_J65A_BHZ_fap01_to_25_calibrated_01_02_10_20. sac 
SAC>  INICM 

Figure  21 .  SAC  Example  of  T ransfer  Command  Using  FAP  File  to 

Calibrate  Data. 


4.  RSAC.m 

The  rsac.m  program  is  a  publicly  available  MATLAB  program  that  reads  a 
SAC  file  and  outputs  a  three  column  MATLAB  file  that  is  easily  useable  for 
analysis.  The  first  rsac.m  column  in  the  output  is  time,  the  second  column  is 
amplitude,  and  the  third  column  contains  the  SAC  header  information  and  was 
not  used  in  this  thesis. 

5.  MATLAB  Data  Analysis 

MATLAB  was  used  for  data  analysis  on  the  calibrated  SAC  files  after  the 

data  was  read  using  rsac.m.  Three  functions  (see  Appendix  A)  were  written  that 

take  as  inputs  a  maximum  desired  amplitude,  the  sampling  frequency  of  the 

OBS,  and  the  data  from  the  three  orthogonal  seismometer  channels  from  an 

OBS.  The  first  function  entitled  “vibration_analysis.m”  finds  the  peaks  of  the 
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signal  above  the  maximum  amplitude  and  checks  the  times  between  these 
peaks.  It  then  looks  for  the  longest  period  of  time  between  peaks.  If  a  quiescent 
period  (a  period  of  time  during  which  the  signal  remains  below  the  maximum 
amplitude)  is  found  in  the  vertical  channel  that  is  at  least  an  hour  long,  a  power 
spectral  density  (PSD)  using  the  pwelch  function  in  MATLAB  is  conducted  on  the 
quiescent  period.  The  PSD  displays  all  three  orthogonal  channels.  The  data 
within  the  taper  at  the  beginning  and  end  of  the  input  time  series  is  not  included 
in  the  quiescent  period  and  is  removed  prior  to  a  PSD  being  conducted.  If  there 
is  not  at  least  an  hour  of  quiescent  time,  a  different  MATLAB  function  was  used. 
If  the  signal  was  too  high  and  no  hour-long  quiescent  periods  exist,  the  function 
high_amp_vibration_analysis.m  was  used.  This  program  finds  the  three-hour 
time  period  with  the  lowest  mean  amplitude  and  outputs  the  same  information  as 
the  vibration_analysis.m  program.  If  the  24-hour  sample  was  low  amplitude  and 
did  not  cross  the  maximum  amplitude  line,  the  low_amp_vibration_analysis.m 
was  conducted,  which  outputs  a  PSD  for  the  entire  sample. 

There  are  also  three  different  types  of  histograms  produced  by  these 
programs.  The  first  is  a  histogram  of  the  amplitudes  throughout  a  24-hour  period. 
Both  non-logarithmic  and  logarithmic  y-axes  were  required  to  adequately  display 
all  of  the  occurrences.  The  next  two  histograms  involved  the  desired  maximum 
amplitude  threshold  which  is  entered  into  the  functions.  The  first  is  a  histogram  of 
width  of  pulses  that  exceed  the  maximum  amplitude  threshold,  and  the  second  is 
a  histogram  of  lengths  of  time  between  these  pulses  that  exceed  the  maximum 
amplitude  threshold.  Both  of  these  histograms  are  completed  logarithmically  and 
non-logarithmically  to  ensure  that  outliers  are  displayed. 

The  maximum  amplitude  threshold  histograms  were  completed  because 
vibration  amplitudes  above  a  certain  level  would  likely  result  in  a  noise  floor  too 
high  to  obtain  satisfactory  performance  for  the  bottom  mounted  sensor.  However, 
since  the  performance  of  a  future  bottom  mounted  sensor  is  unknown,  the 
maximum  amplitude  threshold  was  chosen  somewhat  arbitrarily  to  be  5x1 0"5  m/s 
based  on  observations  of  high  amplitude  activity.  This  threshold  proved  to  be  a 
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good  indicator  of  high  amplitude,  seismic  activity  (non-quiescent),  and  the 
threshold  was  thus  used  to  compare  the  velocities  measured  at  all  of  the  OBS 
locations.  Only  YM  39  had  standard  deviations  for  24-hour  periods  that  regularly 
exceeded  this  threshold,  and  the  other  three  OBSs,  with  the  exception  of  a 
couple  of  days,  remained  well  below. 

E.  CHALLENGES  FACED  WITH  CALIBRATION 

The  largest  challenge  faced  for  this  thesis  was  obtaining  calibrated  data 
before  analysis  could  be  conducted  in  MATLAB.  [8]  describes  the  LDEO  OBS  as 
able  to  effectively  measure  from  .01  to  20  Hz.  The  LDEO  OBS’s  response  curve, 
illustrated  in  Figure  22,  changes  by  2  orders  of  magnitude  from  .01  to  1  Hz.  Figure 
23  exhibits  the  SIO  OBS’s  response  curve  that  has  nearly  a  flat  response  from  .01 
to  20  Hz.  Because  of  the  LDEO  response  curve,  getting  the  calibrated  data 
becomes  more  complicated  than  just  dividing  the  raw  counts  by  the  sensitivity. 


YM.39..BHZ 

2008-  05-  11T21:25:51  to  2009-05- 30T21:20:11 


Figure  22.  LDEO  OBS  Response  Curve  for  YM  39  Channel  Z. 

Source:  [26]. 
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7D.J65A..BHZ 

2011-10-17T22:31:00  to  2012-07-16T07:21:00 


Figure  23.  SIO  OBS  Response  Curve  from  Channel  Z  of  7D  J65A. 

Source:  [27], 

Due  to  lack  of  familiarity  with  SAC  software,  and  because  the  SAC  manual 
does  not  go  into  detail  about  the  process  behind  calibration,  calibration  was 
initially  conducted  with  MATLAB  using  the  CS  file.  This  proved  to  be  more 
complicated  than  necessary,  especially  after  conversation  with  a  seismologist 
who  stated  that  the  SAC  software  was  quite  able  to  complete  the  calibration 
correction  sufficiently.  There  were  four  ways  attempted  to  complete  the 
calibration  that  involved  the  IRIS  SAC  software  and  the  time  series  Webserver 
(TSWS). 

Figure  24  displays  the  TSWS  webpage.  The  TSWS  allows  a  user  to  pick 
the  network,  station,  channel,  date(s),  time,  and  multiple  calibration/filtration 
options.  Depending  on  the  selection,  it  will  immediately  produce  a  plot  or 
downloadable  data  in  a  user-selected  format.  For  rsac.m  to  work  with  the  TSWS 
data,  “SAC  binary  little-endian”  must  be  selected.  With  the  correct  settings,  the 
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TSWS  provided  a  convenient  way  to  double  check  time  plots.  Without  correct 
settings,  the  TSWS  produced  some  odd  results  that  could  be  misleading. 

To  illustrate,  Figure  25  is  a  TSWS  plot  of  a  full  day  from  YM  39,  with 
frequency  limits  of  .01 -.02-1 0-20  Hz,  the  instrument  correction  applied,  and 
“remove  mean”  selected.  Figure  26  is  a  MATLAB  plot  produced  with  FAP 
calibrated  data  of  the  same  time  period  and  settings.  Both  plots  appear  similar. 


Web  Services  IRISWS  Timeseries  Docs  v.  1  Builder 


URL  Builder  timeseries  v.1 


Service  interface  URL  Builder  Help  Revisions 

Use  this  form  to  build  a  URL  to  the  timeseries  web  service.  Notice  that  as  you  edit  the  form,  the  link  is  automaticaly  updated. 


O  Usage 


Network: 

7D 

Remove  mean: 

O 

Station: 

J65A 

Low-Pass  Filter: 

□  1.0 

Location: 

High-Pass  Filter 

Channel: 

BHZ 

Band-Pass  Filter: 

Start  Time: 

2011-11-15T 12:00:00  fi 

Differentiate: 

End  Time: 

2011-11-15T  12:01:21  fi 

Integrate: 

Correction: 

Perform  Instrument  Correction  ; 

Envelope: 

Frequency 

O  .0V.02-10-20 

Taper: 

Limits: 

Decimate  (samples  per 

Units: 

Scale: 

sec): 

Output: 

SAC  binary  little-endian  ; 

Click  the  link: 

http://service.iris.edu/iriswsAimeseries/1/query?net-7D&sta-J65Mcha-BHZ&start-2011-11-15T12:00:00&end-2011-11-15T12:01:21&freqlimits-.01-.02-10-20& 
demean- true&output”sacbl&k>c”-&correct-true 


This  form  primarily  serves  to  illustrate  how  to  create  ws-timeseries  URLs  and  includes  a  fixed  sequence  of  signal  processing  options.  The  ws- timeseries 
service  performs  all  signal  processing  options  in  the  order  they  appear  within  the  URL.  this  flexibility  is  thus  not  offered  by  the  form;  advanced  users 
may  wish  to  manually  modify  or  create  their  own  request  URLs  with  the  processing  options  in  their  preferred  order. 


Figure  24.  TSWS  Webpage.  Source:  [28]. 
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YM_39J?_BHZ,  3456000  samples,  40  0  SpS,  2008-06-13TOO:00:00.002 


Figure  25.  TSWS  Plot  of  YM  39  on  June  13,  2008,  24  Hours, 
Frequency  Limits  of  .01 -.02-1 0-20  Hz.  Source:  [28]. 


YM  3913  JUN  08  All  Day 


Figure  26.  MATLAB  Plot  of  YM  39  on  June  1 3,  2008,  24  Hours, 
Frequency  Limits  of  .01 -.02-1 0-20  Hz,  FAP  Calibrated. 


36 


However,  if  the  frequency  range  is  changed  slightly  to  0-. 02-1 0-20  Hz,  and 
all  other  settings  remain  unchanged,  the  TSWS  produces  the  plot  in  Figure  27. 
Figure  28  displays  the  FAP  calibrated  file  with  the  same  settings.  Although  the 
plot  has  changed  with  apparent  ringing,  it  still  is  recognizable  as  a  seismogram 
with  close  to  the  same  values  as  the  previous  FAP  plot. 


YM.39.7P.8HZ,  3456000  samples,  40.0  sps,  2008-06-13TOO:00:00.002 


Figure  27.  TSWS  Plot  of  YM  39  on  June  13,  2008,  24  Hours, 
Frequency  Limits  of  0-. 02-1 0-20  Hz.  Source:  [28]. 


37 


YM  3913  JUN  08  All  Day 


Figure  28.  MATLAB  Plot  of  YM  39  on  June  1 3,  2008,  24  Hours, 

Frequency  Limits  of  0-. 02-1 0-20  Hz,  FAP  Calibrated. 

Smaller  samples  of  the  same  day  were  inspected  with  similar  results. 
Figure  29  shows  a  TSWS  output  of  4080  samples,  from  12:00:00  to  12:01 :42  on 
13  June  2008.  The  settings  were  the  same  as  for  Figure  25,  with  frequency  limits 
of  .01 -.02-1 0-20  Hz,  instrument  correction  applied,  and  “remove  mean”  selected. 
Figure  30  displays  what  happens  to  this  smaller  data  sample  when  “remove 
mean”  is  not  selected  with  all  other  settings  unchanged. 
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(x  10-5) 


YM_39_??_BHZ,  4080  samples,  40  0  sps,  2008-06-13T12  0000  001 


Figure  29.  TSWS  Plot  of  YM  39  on  June  13,  2008,  from  12:00:00  to 
1 2:01 :42,  Frequency  Limits  of  .01  -.02-1 0-20  Hz.  Source:  [28]. 


YM.39.7P.BHZ,  4080  samples,  40.0  sps.  2008-06-13T12  00  00  001 


Figure  30.  TSWS  Plot  of  YM  39  on  June  13,  2008,  from  12:00:00  to 
12:01 :42,  Frequency  Limits  of  .01 -.02-1 0-20  Hz,  Remove  Mean 
Deselected.  Source:  [28]. 
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Figure  31  is  a  FAP  calibrated  plot  of  the  same  time  period  without  “remove 
trend”  applied.  This  behavior  of  the  TSWS  led  to  the  conclusion  that  other 
calibration  methods  were  needed  besides  using  the  user-friendly  TSWS.  If  the 
frequency  limits  are  applied  correctly,  the  instrument  correction  is  applied,  and 
“remove  mean”  is  selected,  the  amplitude  of  TSWS  plots  can  still  differ  from  the 
FAP,  PZ,  and  RESP  calibrated  plots,  but  the  shape  of  the  signal  is  nearly 
identical  to  the  other  calibration  methods. 


YM  39  13  JUN  08  1200:00-1201:42 


Figure  31 .  MATLAB  Plot  of  YM  39  on  June  1 3,  2008,  1 2:00:00  to 
1 2:01 :42,  Frequency  Limits  of  .01  -.02-1 0-20  Hz,  FAP 
Calibrated,  Remove  Trend  Not  Completed. 


The  remaining  SAC  calibration  methods  produced  very  similar  results  if 
the  frequency  limits  were  set  the  same.  It  was  initially  thought  that  FAP  produced 
superior  results  to  both  a  RESP  and  PZ  calibrated  file,  but  it  was  later  discovered 
that  FAP,  RESP,  and  PZ  calibrated  files  are  nearly  identical.  Figures  32  and  33 
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show  four  different  methods  to  gain  calibrated  data  for  both  OBSs.  Figure  32  is 
from  YM  39  on  June  13,  2008,  and  Figure  33  is  from  7D  J65A  on  November  15, 
2011.  Of  note,  all  the  plots  are  nearly  identical,  but  the  TSWS  produces  different 
amplitudes  in  Figure  32. 


YM  39  13  JUN  08  12:00:00-12:01:42 
RESP  Calibrated 


YM  39  13  JUN  08  1 2:00:00-1 201:42 


YM  39  13  JUN  08  12:00:00-12:01:42 


YM  39  13  JUN  08  12:00:00-12:01:42 


Figure  32.  MATLAB  Plot  of  Different  Calibration  Methods  for 

YM  39  on  June  13,  2008. 
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7D  J65A  15  NOV  2011  1 2:00:00-1 2:01:21 


Time  (s)  x  to4 


7D  J65A  1 5  NOV  201 1  1 2:00:00-1 2:01 :21 


7D  J65A15  NOV  2011  1 2:00:00-1 2:01:21  7D  J65A  1 5  NOV  201 1  1 2:00:00- 1201:21 


Figure  33.  MATLAB  Plot  of  Different  Calibration  Methods  for 
7D  J65A  on  November  1 5,  201 1 . 

Because  of  the  similarity  of  the  data  from  the  three  SAC  calibration 
techniques  and  because  of  the  work  already  completed  with  FAP  files,  the  FAP 
calibration  technique  was  used  for  the  data  in  this  thesis.  Any  of  the  methods 
available  with  the  SAC  transfer,  either  FAP,  RESP,  or  PZ  files,  would  have  been 
sufficient  for  use. 

Some  idiosyncrasies  with  the  software  proved  to  hold  additional 
challenges.  Each  time  a  transfer  command  in  SAC  was  used,  the  command  line 
involved  a  “from”  type  of  file  (i.e.  RESP,  PZ,  FAP)  and  an  optional  “to”  desired 
units  (e.g.,  velocity  and  displacement).  Figure  21  is  an  example  of  transfer  with  a 
“to  vel”  (to  velocity)  command.  The  default  waveform  is  displacement,  so  if  a  “to 
none”  line  is  entered,  the  SAC  manual  [24]  states  that  the  output  will  default  to 
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displacement.  If  no  “to”  line  is  entered  (left  blank),  SAC  will  assume  it  is  “to  none” 
[24],  A  transfer  with  a  FAP  file,  though,  does  not  necessarily  follow  the  SAC 
manual,  which  led  to  an  assumption  that  calibrated  files  were  in  units  of  velocity 
when  the  files  output  were  actually  acceleration.  Figure  34  summarizes  the  SAC 
peculiarities  with  the  transfer  command  that  led  to  some  analysis  challenges. 


POLEZERO  VEL  DISP  DISP 

FAP  ACC  VEL  VEL 


Figure  34.  Summary  of  SAC  Transfer  Output. 

F.  SUMMARIZED  DATA  ANALYSIS  STEPS 

The  steps  described  in  this  chapter  are  summarized  in  Figure  35. 


Figure  35.  Methodology  and  Data  Analysis  Summary. 
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V.  RESULTS 


A.  POWER  SPECTRAL  DENSITY 

1.  YM  39  and  YM  38 

YM  39  was  the  most  erratic  of  the  OBSs  sampled.  It  had  the  highest 
variation  in  standard  deviation  and  the  largest  variation  in  percentage  of  time 
spent  above  max  amplitude.  Most  samples  did  not  have  at  least  one  hour  of 
quiescent  time  and  had  to  be  analyzed  with  the  high  amplitude  vibration  analysis 
code.  In  most  instances  for  YM  39,  the  PSDs  taken  from  the  quiescent  periods 
displayed  a  decreasing  amplitude  with  increasing  frequency  and  no  apparent 
narrowband  frequencies  as  illustrated  in  Figure  36;  however,  the  decibel  levels 
differed  from  sample  to  sample.  Discrete  frequencies  do  show  up  in  the  PSDs  of 
the  entire  24-hour  period,  which  is  displayed  in  Figure  37. 


Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude 


Figure  36.  YM  39,  PSD  from  August  1 , 2008,  Quiescent  Period  (Three 

Hours  with  Lowest  Mean  Amplitude). 
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Power  Spectral  Density  Estimate 
Entire  Day  of  Data 


Figure  37.  YM  39,  August  1 , 2008,  PSD  for  Entire  Day 

If  the  PSD  for  channel  Z  is  graphed  on  a  logarithmic  x  and  y  axis,  it  appears 
almost  as  a  straight  line  as  in  Figure  38.  Most  spectral  densities  for  YM  39 
channel  Z  follow  this  pattern  with  only  slight  variations  from  the  overall  slope.  The 
velocity  time  series  plot  responsible  for  this  PSD,  Figure  39,  represents  a  typical 
appearance  of  the  majority  of  the  days  for  YM  39. 
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Power  Spectral  Density  Estimate 


Figure  38.  YM  39,  August  1 , 2008,  PSD  with  Logarithmic  Axes. 


Figure  39.  YM  39  Velocity  Time  Series,  August  1 , 2008. 

For  the  low  frequency  range,  shown  in  Figure  40,  most  of  YM  39’s  PSDs 
exhibited  a  steady  decline  with  the  same  overall  pattern  from  month  to  month. 
However,  even  though  the  overall  shape  and  slope  remained  consistent,  the 
decibel  level  varied  greatly  from  sample  to  sample  and  was  far  more  inconsistent 
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than  the  other  OBSs.  The  drop  off  in  decibel  levels  close  to  0  Hz  is  due  to  the 
bandpass  filter,  which  starts  with  a  taper  at  .01  Hz. 


Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude,  0  to  2  Hz 


Figure  40.  YM  39,  PSD,  August  1 , 2008,  0  to  2  Hz. 

YM  38,  however,  had  the  lowest  amplitudes  of  all  four  OBSs,  especially  in 
the  vertical  channel.  Channel  1  and  channel  2  exhibited  much  higher  amplitudes 
than  channel  Z  for  all  days  sampled.  This  is  shown  by  the  fact  that  the  decibel 
level  of  the  vertical  channel  hovered  around  10  to  20  dB  less  than  the  horizontal 
channels.  The  channel  Z  plot  showed  hints  of  discrete  frequencies,  especially 
around  4  Hz,  6  Hz,  and  8  Hz,  but  the  appearance  of  the  discrete  frequencies 
differed  slightly  from  month  to  month  like  in  Figure  41  where  the  discrete 
frequencies  are  located  at  approximately  4  Hz  and  7  Hz.  A  visual  comparison  of 
plots  between  YM  38  and  YM  39  on  the  same  dates  revealed  no  common 
discrete  frequencies  and  different  appearances  of  the  time  series  plots.  In 
particular,  YM  38  was  quieter  than  YM  39  by  around  2  orders  of  magnitude. 
Figure  41  displays  a  typical  YM  38  PSD  for  a  quiescent  period.  Figure  42  is  from 
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the  same  sample  between  0  to  2  Hz.  Most  24-hour  samples  from  YM  38  never 
had  any  amplitude  spikes  above  the  max  amplitude  threshold. 


Power  Spectral  Density  Estimate 
Longest  Continuous  Quiescent  Period 


Figure  41 .  YM  38,  PSD  for  Quiescent  Period,  August  1 , 2008. 


Power  Spectral  Density  Estimate 
Longest  Continuous  Quiescent  Period,  0  to  2  Hz 


Frequency  (Hz) 


Figure  42.  YM  38,  PSD  for  Quiescent  Period,  August  1 , 2008,  0  to  2 

Hz. 
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The  sample  from  August  1,  2008  is  an  example  of  one  of  the  few  days  when 
the  amplitude  crossed  the  maximum  amplitude  threshold.  As  displayed  in  Figure 
43,  there  are  two  large  amplitude  occurrences.  Figure  44  is  a  close-up  of  the 
largest  amplitude  spike,  and  the  amplitude  increase  does  not  appear  to  cause 
any  ringing  on  the  plot. 


Figure  43.  YM  38,  Velocity  Time  Series,  August  1 , 2008. 


Close-up  of  High  Amplitude.  1  AUG  2008,  YM  38.  Vertical  Velocity  (Channel  Z) 


Figure  44.  YM  38,  August  1 , 2008,  Close-Up  of  Highest  Amplitude 

Event  in  Velocity  Time  Series 
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Unlike  in  OBS  YM  39,  most  of  the  high  amplitudes  in  the  horizontal 
channels  do  not  show  up  in  the  vertical  channel.  Figure  45  is  a  depiction  of 
channel  1  from  YM  38  on  August  1 , 2008.  The  higher  amplitudes  in  the  beginning 
of  the  day  do  not  translate  into  a  high  response  in  channel  Z,  but  the  two  high- 
amplitude  occurrences  after  the  3x1 04  second  mark  and  4x1 04  second  mark  do 
appear  in  the  vertical  channel. 


Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity  (Channel  1) 


Tim#  M 


Figure  45.  YM  38  Velocity  Time  Series,  Channel  1 ,  August  1 , 2008. 

Although  YM  38  and  YM  39  were  located  in  the  Luzon  Strait  in  close 
proximity,  the  two  OBSs  displayed  vastly  different  results.  YM  39  was  not 
consistent,  and  YM  38  displayed  much  lower  vertical  vibrational  disturbances. 

2.  7D  J65A  and  7D  J73A 

Both  OBSs  west  of  the  Strait  of  Juan  de  Fuca,  7D  J65A  and  7D  J73A 
display  similar  and  consistent  results.  7D  J65A  exhibits  the  same  consistent 
pattern  almost  every  day  where  channel  Z  peaked  at  8  to  9  Hz  with  varying 
decibel  levels.  Figures  46  and  47  depict  the  two  7D  J65A  plots  with  maximum 
and  minimum  decibel  levels  at  around  9  Hz.  Channel  Z  in  Figure  46  peaks  at 
around  -130  dB  and  in  Figure  47  peaks  at  around  -150  dB. 
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Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude 


Figure  46.  7D  J65A,  PSD  of  Quiescent  Period  (Three  Hours  with 
Lowest  Mean  Amplitude),  February  15,  2012. 


Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude 


Figure  47.  7D  J65A,  PSD  of  Quiescent  Period  (Three  Hours  with 
Lowest  Mean  Amplitude),  March  1, 2012. 
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Figure  48  is  from  7D  J73A,  March  1,  2012,  and  shows  an  increase  in 
decibels  after  10  Hz.  Even  though  a  low-pass  taper  is  applied  at  10  Hz,  7D 
J73A’s  decibel  levels  continue  to  increase  after  this  taper  is  applied.  Decibel 
levels  for  channel  Z  peak  after  the  taper  from  as  low  as  -160  dB  to  as  high 
as  -130  dB.  This  pattern  of  decibel  level  increase  after  the  taper  is  present  in 
nearly  every  sample. 


Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude 


Figure  48.  7D  J73A,  PSD  of  Quiescent  Period  (Three  Hours  with 
Lowest  Mean  Amplitude),  March  1, 2012. 


From  0  to  2  Hz,  7D  J65A  and  7D  J73A  are  similar,  with  both  OBSs  peaking 
at  around  0.2  Hz  and  then  displaying  a  steady,  almost  linear,  decline  after  the 
peak.  Figures  49  and  50  depict  0  to  2  Hz  from  both  OBSs  on  the  same  day, 
March  1, 2012. 
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Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude,  0  to  2  Hz 


Figure  49.  7D  J65A,  PSD  of  Quiescent  Period  (Three  Hours  with 
Lowest  Mean  Amplitude),  March  1 , 201 2,  0  to  2  Hz. 


Power  Spectral  Density  Estimate 
Three  Hours  with  Lowest  Mean  Amplitude,  0  to  2  Hz 


Figure  50.  7D  J73A,  PSD  of  Quiescent  Period  (Three  Hours  with 
Lowest  Mean  Amplitude),  March  1 , 201 2,  0  to  2  Hz. 
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3. 


PSD  Conclusion 


Figures  51  and  52  summarize  the  OBS  PSD  variability,  focusing  in  on  1 
Hz  and  5  Hz.  The  decibel  levels  for  both  1  Hz  and  5  Hz  are  graphed  per  sample 
for  all  samples  taken.  The  results  show  fairly  steady  decibel  levels  for  all  OBSs  at 
the  frequencies  except  YM  39,  which  is  the  least  consistent  and  noisiest  of  all  the 
OBSs  sampled.  The  decibel  levels  at  5  Hz  are  similar  for  YM  38  and  the  two 
OBSs  located  west  of  the  Strait  of  Juan  de  Fuca.  At  1  Hz,  though,  7D  J65A  and 
7D  J73A  are  consistently  stronger.  Tables  listing  all  the  values  in  Figures  51  and 
52  are  located  in  Annex  B. 


PSD  Values  at  1  Hz  and  5  Hz,  YM  39  and  YM  38 
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Figure  51 .  Quiescent  PSD  Values  at  1  Hz  and  5  Hz  for 

YM  39  and  YM  38. 
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Figure  52.  Quiescent  PSD  Values  at  1  Hz  and  5  Hz  for 

7D  J65A  and  J73A. 


55 


Although  acceleration  is  not  discussed  in  Results,  acceleration  can  be 
obtained  from  the  velocity  data  used  in  all  of  the  figures  by  using  Equation  21, 
where  a  is  acceleration,  co  is  angular  frequency,  and  v  is  velocity.  Appendix  B 
contains  rms  acceleration  for  all  vertical  channels. 

a  =  icov  (21) 

B.  AMPLITUDE  HISTOGRAMS 
1.  7D  J65A  and  7D  J73A 

Amplitude  histograms  were  generally  viewed  with  a  logarithmic  y-axis  so 
that  more  details  were  available.  However,  any  velocity  time  series  plot  that  was 
consistent  and  without  multiple  large  amplitude  spikes  appeared  Gaussian  if 
viewed  without  a  logarithmic  axis.  Most  of  the  7D  J65A  and  7D  J73A  amplitude 
histograms  were  Gaussian  in  appearance.  Every  amplitude  histogram  from  these 
two  OBSs  viewed  on  a  logarithmic  scale  looked  similar  to  either  Figure  53  or 
Figure  56.  Figure  53  is  from  OBS  7D  J73A,  November  1,  2011.  Figure  54 
displays  the  same  histogram,  but  with  a  non-logarithmic  y-axis.  The  channel  Z 
time  series  plot  represented  by  these  histograms  is  shown  in  Figure  55. 
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Histogram  of  Channel  Z  Amplitudes 


Figure  53.  7D  J73A,  Amplitude  Histogram,  November  1 , 201 1 . 


-1.5  -1  -0.5  0  05  1  15  2 

Velocity  Amplitudes  (nVs)  xio-4 


Figure  54.  7D  J73A,  Amplitude  Histogram,  November  1 , 201 1 ,  Non- 

Logarithmic. 
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Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Vertical  Velocity  (Channel  Z) 
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Figure  55.  7D  J73A,  Velocity  Time  Series,  November  1 , 201 1 . 

Figure  56  displays  a  narrower  histogram  in  a  logarithmic  plot.  Such  narrow 
histograms  were  also  common  in  many  of  the  days  from  7D  J65A  and  7D  J73A. 
The  histogram  in  Figure  56  naturally  has  a  smaller  number  of  amplitude  events 
that  surpass  the  amplitude  threshold.  Figures  57  and  58  display  the  same 
information  in  non-logarithmic  plots,  but  Figure  58  is  a  close-up  displaying  the 
single,  larger  amplitude  events  that  occurred  throughout  this  sample  day.  Figure 
59  is  the  velocity  time  series  plot  of  channel  Z  for  July  1,  2012,  represented  by 
the  histograms  in  Figures  56  through  58.  Both  OBSs  are  consistent  enough  that 
the  wider  and  narrow  histograms  in  Figure  53  and  Figure  56  are  representative  of 
the  range  of  standard  deviation  and  shapes  of  all  the  data  investigated. 


58 


Figure  56.  7D  J73A,  Amplitude  Histogram,  July  1 , 2012. 
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Figure  57.  7D  J73A,  Amplitude  Histogram,  July  1 , 201 2,  Non- 

Logarithmic. 
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Figure  58.  7D  J73A,  Amplitude  Histogram,  July  1 , 2012, 

Close-Up  Non-Logarithmic. 
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Figure  59.  7D  J73A,  Velocity  Time  Series,  July  1 , 201 2, 

Channel  Z. 


2.  YM  39  and  YM  38 

For  the  YM  39  and  YM  38  OBSs,  the  non-logarithmic  amplitude 
histograms  miss  many  details  otherwise  available  with  a  logarithmic  plot.  Many 
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time  series  plots  for  both  YM  38  and  YM  39  involve  long,  relatively  quiet  time 
periods  followed  by  shorter  periods  of  higher  amplitude  vibrations.  Figure  60  from 
June  13,  2008,  displays  this  type  of  time  series  plot.  The  amplitude  histogram  in 
Figure  61  is  logarithmic,  and  the  amplitude  histogram  in  Figure  62  is  not  and 
misses  useful  specifics  provided  in  Figure  61 .  Many  of  the  histograms  for  YM  38 
and  YM  39  have  similar  appearances,  although  the  samples  from  YM  38  were  up 
to  two  orders  of  magnitude  smaller  in  standard  deviation. 


Figure  60.  YM39,  Velocity  Time  Series,  June  13,  2008,  Channel  Z. 


61 


Figure  61 .  YM  39,  Amplitude  Histogram,  June  1 3,  2008. 


Figure  62.  YM  39  Amplitude  Histogram,  Non-Logarithmic, 

June  13,2008. 
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Multiple  amplitude  histograms  from  both  YM  38  and  YM  39  appear  in  a 
crown-type  shape  shown  in  Figure  63.  Figure  64  is  the  time  series  plot  for  the 
same  day,  February  1,  2009.  This  histogram  illustrates  the  common  shape  for 
the  plots  when  there  were  calm  periods  interrupted  by  multiple,  almost  impulsive, 
amplitude  spikes.  Figures  65  and  66,  from  August  1,  2008,  display  the  same 
behavior  with  YM  38  but  at  a  much  smaller  amplitude.  Initially  this  behavior  was 
thought  to  be  electronic  clipping,  but  there  are  histogram  examples  with  the  same 
appearance  but  increased  amplitude.  The  vibrations  arrive  and  are  symmetric 
with  an  equal  positive  and  negative  amplitude  response. 


Histogram  of  Channel  Z  Amplitudes 
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Figure  63.  YM  39,  Amplitude  Histogram,  February  1, 2009. 
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Figure  64.  YM  39,  Velocity  Time  Series,  February  1 , 2009. 
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Figure  65.  YM  38,  Amplitude  Histogram,  August  1 , 2008. 
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Figure  66.  YM  38,  Velocity  Time  Series,  August  1 , 2008,  Channel  Z. 

3.  Amplitude  Histogram  Conclusion 

The  histograms  for  both  the  OBSs  located  west  of  the  Strait  of  Juan  de 
Fuca  proved  that  the  data  these  two  OBSs  collected  was  consistent  with  few 
outliers.  This  led  to  amplitude  histograms  with  the  appearance  of  a  normal 
distribution.  Many  samples  taken  from  the  OBSs  in  the  Luzon  Strait  displayed 
long,  relatively  quiet  periods,  which  were  abruptly  interrupted  by  higher  amplitude 
seismic  readings.  These  particular  events  were  responsible  for  the  crown-type 
shapes  in  Figures  63  and  65. 

C.  HISTOGRAMS  OF  HIGH  AMPLITUDE  PULSE  WIDTHS  AND 

HISTOGRAMS  OF  TIME  BETWEEN  HIGH  AMPLITUDE  PULSES 

Two  types  of  histograms  were  based  on  a  maximum  entered  amplitude 
threshold:  a  histogram  of  the  width  of  pulses  that  exceeded  the  maximum 
amplitude  and  a  histogram  of  times  between  these  pulses.  As  mentioned  in 
Chapter  IV,  the  maximum  amplitude  threshold  chosen  for  this  thesis  was  5x1 0"5 
m/s.  Although  this  threshold  was  higher  than  the  standard  deviation  of  any  OBS 
except  for  YM  39,  all  of  the  OBSs  have  occurrences  where  the  amplitude 
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surpasses  this  line.  This  threshold  proved  to  be  a  good  indicator  of  high 
amplitude,  non-quiescent  activity. 

1.  7D  J65A  and  7D  J73A 

Both  OBSs  west  of  the  Strait  of  Juan  de  Fuca  spent  around  6%  of  the 
sampled  time  above  the  maximum  amplitude  threshold.  For  these  OBSs,  the 
maximum  signal  length  above  the  maximum  amplitude  threshold  was  around  6 
seconds,  but  the  majority  of  these  high  amplitude  signals  were  1  second  or  less 
in  length.  Figure  67  is  a  representative  histogram  of  signals  that  exceeded  the 
maximum  amplitude  threshold  for  both  7D  J65A  and  7D  J73A.  Figure  68  is  the 
time  series  associated  with  this  histogram. 


Histogram  of  Length  of  Continuous  Time  Signal  Exceeding  Maximum 
Channel  Z 
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Figure  67.  7D  J65A,  November  1 5,  201 1 ,  Width  of  Signals 
Exceeding  Maximum  Amplitude  Threshold 
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Figure  68.  7D  J65A,  November  1 5,  201 1 , 

Velocity  Time  Series 

The  vast  majority  of  times  between  signals  that  exceeded  the  maximum 
amplitude  was  1  second  or  less.  Figures  69  and  70  display  time  between  events, 
and  Figure  71  is  a  close  up  with  a  bin  width  of  1  second.  The  higher  amplitude 
signals  occur  in  groups,  which  explains  the  large  number  of  occurrences  of  the  1 
second  or  less  bin. 


Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude 
Channel  Z 
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Figure  69.  7D  J65A,  November  1 5,  201 1 ,  Time  Between 
Signals  That  Exceed  Maximum  Amplitude 
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Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude 
Channel  Z 
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Figure  70.  7D  J65A,  November  1 5,  201 1 ,  Time  Between  Signals  That 
Exceed  Maximum  Amplitude,  Non-Logarithmic 
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Figure  71 .  7D  J65A,  November  1 5,  201 1 ,  Time  Between 
Signals  That  Exceed  Maximum  Amplitude 
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2. 


YM  39  AND  YM  38 


OBS  YM  39  spent  an  average  of  52%  of  time  above  the  maximum 
amplitude  threshold  over  all  the  days  sampled,  although  the  percentage  of  time 
spent  above  the  threshold  varied  significantly  from  1%  to  91%  per  day  sampled. 
The  sample  from  March  15,  2009,  spent  54%  of  the  time  above  the  threshold, 
which  is  close  to  the  overall  YM  39  average  of  52%.  Figure  72  depicts  the  time 
series  of  channel  Z  from  March  15,  2009.  Figure  73  is  the  histogram  of  signal 
widths  above  the  maximum  amplitude  threshold,  and  Figure  74  displays  the 
times  between  these  signals.  Figure  75  is  non-logarithmic,  and  displays  the 
single  occurrences  where  the  largest  time  period  between  signals  is  almost  five 
minutes. 


Figure  72.  YM  39,  March  15,  2009,  Velocity  Time  Series 
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Histogram  of  Length  of  Continuous  Time  Signal  Exceeding  Maximum 
Channel  Z 
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Figure  73.  YM  39,  March  15,  2009,  Width  of  Signals  Exceeding 

Maximum  Amplitude  Threshold 


Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude 
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Figure  74.  YM  39,  March  15,  2009,  Time  Between  Signals  That  Exceed 

Maximum  Amplitude  Threshold 
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Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude 

Channel  Z 
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Figure  75.  YM  39,  March  15,  2009,  Time  Between  Signals 
That  Exceed  Maximum  Amplitude  Threshold 

3.  Conclusion 

It  was  anticipated  that  the  length  of  time  between  high  amplitude  signals 
would  be  short  for  the  days  with  a  high  percentage  of  time  spent  above  the 
maximum  amplitude  threshold.  What  was  unexpected,  however,  was  the  fact  that 
the  majority  of  times  between  high  amplitude  signals  were  1  second  or  less  for 
an  OBS  that  spent  only  7%  of  its  time  above  the  maximum  amplitude  (Figure  68). 
This  confirms  that  the  high  amplitude  signals  tend  to  occur  in  groups.  The 
maximum  time  spent  between  signals  on  this  same  day  was  16  minutes,  but  this 
only  occurred  once  in  a  24  hour  sample.  Table  2  summarizes  the  extremes  of  the 
OBSs. 
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Table  2.  Maximum  and  Minimum  Standard  Deviations  and  Times 
Spent  Above  the  Maximum  Amplitude  for  All  OBSs 


OBS 

Max  standard  deviation 
observed  (pm/s) 

Min  standard  deviation 
observed  (pm/s) 

Max  %  of  time  spent 
above  maxamplitude 

Min  %  of  time  spent 
above  max  amplitude 

YM  39 

662 

16.1 

91 

1 

YM  38 

7.3 

0.23 

0.16 

0 

7DJ65A 

117 

7.2 

24 

0 

7DJ73A 

39.9 

7.6 

21 

0 
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VI.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  METHODOLOGY 

The  global  view  on  the  IRIS  website  proved  useful  in  the  selection  of 
OBSs  that  were  co-located  in  significant  geographic  chokepoints.  IRIS  also  made 
it  convenient  to  select  the  OBS  directly,  view  the  detailed  information,  and 
request  the  data  directly  from  the  OBS  information  webpage.  The  publicly 
available  seismic  processing  software  on  the  IRIS  website  was  essential  in 
obtaining  calibrated  data  from  the  raw  data.  Manual  calibration  was  attempted, 
but  the  IRIS  software,  which  was  designed  by  the  seismology  community,  was 
quite  capable  once  familiarity  with  the  software  was  established.  There  were 
multiple  paths  to  calibrated  data,  especially  using  the  SAC  transfer  function,  but 
the  novice  user  will  need  to  ensure  that  the  SAC  transfer  output  is  in  the  units 
that  he/she  is  expecting.  A  Mac  computer,  or  a  Unix-based  operating  system,  is 
required  to  operate  the  IRIS  software.  The  MATLAB  rsac.m  function  was  a 
convenient  method  to  convert  SAC  data  to  a  MATLAB  format. 

B.  RESULTS 

Out  of  the  four  OBSs,  three  proved  consistent  from  month  to  month  in 
their  individual  discrete  frequencies  and  decibel  levels,  but  there  were  no 
common  discrete  frequencies  between  OBSs  that  were  located  in  the  same 
geographic  area.  YM  39  in  the  Luzon  Strait  was  erratic  with  large  differences  in 
decibel  levels  from  sample  to  sample.  The  amplitude  histograms  for  7D  J65A  and 
7D  J73A  displayed  what  appeared  to  be  normal  distributions,  with  the  exception 
of  a  couple  of  outliers  with  some  seismic  activity  noted.  YM  39  and  YM  38 
displayed  oddly  shaped  amplitude  histograms  that  were  the  result  of  long 
quiescent  periods  interrupted  by  large,  impulsive  seismic  activity.  The  histograms 
based  on  the  maximum  amplitude  threshold  showed  that  the  higher  amplitude 
vibrations  occurred  in  groups,  with  the  vast  majority  of  counts  occurring  with  a 
second  or  less  between  events.  Even  the  histograms  from  samples  where  the 
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percentage  of  time  spent  above  the  maximum  amplitude  threshold  was  low 
showed  results  where  the  vast  majority  of  time  between  signals  was  one  second 
or  less. 

The  extreme  differences  between  YM  38  and  YM  39  lead  to  questions 
about  environmental  effects  on  YM  39  that  are  not  affecting  YM  38.  Even  though 
YM  38  exhibited  some  high-amplitude  activity  on  occasion,  it  was  far  more  quiet 
and  consistent  than  YM  39.  There  was  not  much  in  the  way  of  commonality 
between  these  two  OBSs  despite  their  close  proximity  in  the  same  geographic 
area. 

7D  J65A  and  7D  J73A  had  close  values  in  decibel  levels,  standard 
deviation,  and  percentage  of  time  spent  above  the  maximum  amplitude  threshold 
in  almost  every  sample  taken  on  the  same  day,  with  the  exception  of  a  couple  of 
outliers  higher  amplitude  seismic  events.  It  would  appear  that  west  of  the  Strait  of 
Juan  de  Fuca  is  a  fairly  homogeneous  location  with  regard  to  environmental 
factors  on  the  ocean  floor. 

C.  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

1.  Transfer  Functions  to  Reduce  Long  Period  Noise 

Pressure  from  infragravity  waves  on  the  surface  of  the  ocean  can 
penetrate  to  the  ocean  floor  causing  deformation  of  the  ocean  floor  [5].  According 
to  [5],  these  pressure  variations  are  the  source  of  long  period  seismic  noise,  from 
around  30  seconds  to  100  or  more  seconds  (1  to  3.3  mHz).  These  pressure 
variations  are  also  picked  up  by  the  DPG.  Therefore,  a  transfer  function  can  be 
developed  that  can  get  rid  of  the  seismic  noise  caused  by  the  surface  waves. 
This  could  also  possibly  be  used  to  take  noise  out  of  any  future  bottom  mounted 
sensor. 

2.  Environmental  Effects  in  the  Luzon  Strait 

The  extreme  differences  between  the  two  OBSs  in  the  Luzon  Strait 
present  questions  about  environmental  effects  on  YM  39  that  are  not  affecting 
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YM  38.  The  Luzon  Strait  is  known  for  its  large  internal  wave  activity,  and 
investigating  a  correlation  between  this  activity  with  OBS  results  would  be  an 
interesting  study.  A  detailed  study  of  currents  and  tides  in  the  local  region 
combined  with  OBS  results  could  also  prove  useful. 

3.  Correlation  of  Data  with  Other  OBS  in  Luzon  Strait,  Possible 
Mechanical  Failure  of  YM  39 

Other  OBSs  in  the  neighborhood  of  YM  38  and  YM  39  could  be  inspected 
to  correlate  the  results  obtained  in  this  thesis.  IRIS  data  problem  reports  were 
searched,  but  no  OBSs  in  the  YM  network  were  listed.  Currently,  then,  there  are 
no  indications  of  a  mechanical  failure  of  YM  39,  but  additional  data  could  be 
explored  to  both  confirm  the  results  and  rule  out  any  mechanical  failure  of  YM  39. 

4.  Electronic  Noise  Floor  of  OBSs 

The  electronic  noise  floors  of  the  OBSs  were  not  found  during  the 
research  for  this  thesis,  and  it  is  possible  that  the  amplitude  histograms  with  a 
Gaussian  appearance  were  the  result  of  OBS  electronic  noise.  Thus,  the  OBSs 
could  be  limited  by  the  electronic  noise  floor  during  quiescent  periods.  No 
instrument  noise  models  were  found  that  provide  the  rough  order  of  magnitude  of 
the  electronic  noise  floor,  and  any  future  study  should  explore  these  questions 
further.  In  particular,  it  would  be  valuable  to  have  confidence  in  the  ability  of  OBS 
sensors  to  accurately  measure  ocean  accelerations  down  to  10"6  m/s2  without 
running  into  the  electronic  noise  floor. 

5.  OBS  Requests 

The  OBSIP  website  [4]  contains  information  for  researchers  to  request 
OBSIP  instrument  usage,  and  the  OBSIP  allows  their  OBS  equipment  to  be 
available  to  other  government  research  or  educational  institutions  [4],  Therefore, 
LLNL  or  NPS  could  request  usage  of  OBSs  in  other  areas  of  interest. 
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APPENDIX  A.  MATLAB  CODE 


A.  VIBRATION  ANALYSIS  CODE 

function  vibration_analysis(rsac_channel_Z,max_amp,Fs,rsac_channel_1  ,rsac_channel_2) 
%  OBS,  ocean  bottom  seismometer,  maximum  amplitude,  BHZ,  sampling  frequency 

% 

%  Purpose:  Takes  three  calibrated  channels  of  OBS  data,  OBS  sampling 
%  frequency,  and  maximum  desired  amplitude,  and  finds  the  maximum 
%  quiescent  period  (at  least  1  hour  in  length  where  the  amplitude 
%  does  not  exceed  maximum  amplitude).  Function  outputs  an  amplitude  (m/s) 

%  vs  time  plot  (for  all  three  channels),  the  mean  and  standard  deviation 
%  of  entire  sample  amplitude,  power  spectral  densities  of  longest 
%  quiescent  period,  a  percentage  of  time  spent  above  maximum  desired 
%  amplitude,  a  histogram  of  amplitudes,  a  histogram  of  lengths  of  signals 
%  greater  than  max  desired  amplitude,  and  a  histogram  of  time  periods 
%  between  signals  that  exceed  maximum  amplitude. 

% 

%  Output  Variables: 

%  std_z  =  standard  deviation  of  channel  z  amplitude 
%  percentage_time_above_max_amp  =  percentage  of  time  where  amplitude 
%  is  greater  than  max  amplitude 

% 

%  Input  Variables: 

%  rsac_channel_Z  =  vertical  .sac  file  after  it  has  been  read  through 
%  rsac.m 

%  max_amp  =  desired  maximum  amplitude  (m/s) 

%  Fs  =  sampling  frequency  of  OBS  (Hz) 

%  rsac_channel_1  =  channel  1  .sac  file  after  it  has  been  read  through 
%  rsac.m 

%  rsac_channel_2  =  channel  2  .sac  file  after  it  has  been  read  through 
%  rsac.m 
% 

%  Local  Variables: 

%  cal_data  =  calibrated  vertical  channel  amplitude  (m/s) 

%  time  =  time  (seconds) 

%  PKS  =  value  of  peaks  greater  than  max  amplitudel  (m/s) 

%  LOCS  =  vector  of  location  of  PKS  (integer  location) 

%  loopjndex  =  binary  file  (index)  for  location  of  PKS 

%  loopjndexjogical  =  logical  version  of  loopjndex 
%  cal_data1  =  same  as  cal_data,  but  modified  to  zero-out  peak 
%  locations  (m/s) 

%  PKS2  =  value  of  peaks  in  second  peaks-finding  function  (m/s) 

%  LOCS2  =  vector  of  location  of  PKS2 
%  loopJndexJogical2=logical  version  of  modified  loopjndex 
%  cal_data2  =  same  as  cal  datal ,  but  modified  to  zero-out  additional 
%  peaks  (m/s) 

%  cal_data3  =  same  as  cal_data2,  but  modified  to  contain  longest 
%  quiescent  occurrence  only  (m/s) 

%  largest_quiescentjndex  =  index  locating  the  largest  quiescent 
%  period 

%  channel1_cal_data  =  calibrated  amplitude  of  channel  1  (m/s) 
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%  channel2_cal_data  =  calibrated  amplitude  of  channel  2  (m/s) 

%  time_quiescent  =  time  location  of  longest  quiescent  period  (s) 

%  N  =  FFT  length 
%  overlap  =  pwelch  overlap  value 

%  Pxx  =  channel  Z  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  F,  FI ,  F2  =  Frequency  (Flz) 

%  Pxxl  =  channel  1  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  Pxx2  =  channel  2  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  index  =  index  locating  where  cal_data  is  greater  than  the  maximum 
%  amplitude 

%  above_max_amp_time  =  time  where  cal_data  is  greater  than  maximum 
%  amplitude  (s) 

%  outlier  =  logical  index  location  where  Z-amplitude  is  greater  than 
%  maximum  amplitude 

%  cal_data_update  =  same  as  cal_data,  but  only  contains  amplitudes 
%  greater  than  maximum  amplitude,  (m/s) 

%  W  =  vector  with  lengths  of  signals  that  exceed  max  amplitude 
%  S  =  vector  that  contains  different  lengths  of  time  between  signals 
%  that  exceed  max  amplitude 

% 

%  Functions  Called:  None. 

% 

%  Filename:  vibration_analysis.m 
%  Written  by:  Jeremy  R.  Hankins 

cal_data=rsac_channel_Z(:,2); 

time=rsac_channel_Z(:,1); 

cal_data_qplot=cal_data; 

rsac1_qplot=rsac_channel_1(:,2); 

rsac2_qplot=rsac_channel_2(:,2); 

%  plot  amplitude  (m/s)  vs.  time  (s)  for  all  channels 

figure 

plot(time,cal_data) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],'r-') 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
xlabel('Time  (s)','Fontsize',12) 

title('Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Vertical  Velocity  (Channel 
Z)','Fontsize',16) 

legend('Calibrated  amplitude', '+/-  max  amplitude’) 
figure 

plot(time,rsac_channel_1  (:,2)) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],’r-') 
ylabel('Velocity  (m/s)’,'Fontsize',1 2) 
xlabel('Time  (s)','Fontsize’,12) 
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title('Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity 

(Channel  1)','Fontsize',16) 

legend('Calibrated  amplitude', '+/-  max  amplitude’) 

figure 

plot(time,rsac_channel_2(:,2)) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],'r-') 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
xlabel('Time  (s)','Fontsize’,12) 

title(’Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity 

(Channel  2)','Fontsize',16) 

legend('Calibrated  amplitude’, ’+/-  max  amplitude’) 

rsac_channel_1=rsac_channel_1(:,2); 

rsac_channel_2=rsac_channel_2(:,2); 

%  find  the  peaks  above  the  maximum  amplitude 

[PKS,LOCS]=findpeaks(abs(cal_data),'MinPeakHeight’,max_amp); 

%  tor-loop  to  create  index  locating  where  time  between  peaks  is  at  least  an 
%  hour 

for  i=1  :length(LOCS); 
if  i==1 

if  time(LOCS(i))-time(1  )>=(60*60); 

loop_index(1  :LOCS(i))=0; 
else  loop_index(1  :LOCS(i))=1 ; 
end 

elseif  i>1  &&  iclength(LOCS); 

if  time(LOCS(i))-time(LOCS(i-1  ))>=(60*60) 
loop_index(LOCS(i-1):LOCS(i))=0; 
else  loopJndex(LOCS(i-1):LOCS(i))=1 ; 
end 

elseif  i==length(LOCS); 

if  length(time)-time(LOCS(i))>=(60*60) 
loop_index(LOCS(i):length(time))=0; 
else  loopJndex(LOCS(i):length(time))=1 ; 
end 
end 
end 

loopJndex_logical=logical(loopJndex); 

cal_data1  =cal_data; 

cal_data1  (ioopjndex_logical)=0; 

%  get  rid  of  remaining  peaks  (which  are  located  directly  next  to  previous 
%  peak  locations) 
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[PKS2,LOCS2]=findpeaks(abs(cal_data1),'MinPeakHeight',max_amp); 
for  i=1  :length(LOCS2) 

loop_index(LOCS2(i)-60*Fs:LOCS2(i)+60*Fs)=1 ;  %taking  away  1  minute  on  either  side  of 
remaining  peaks 
end 

loopJndex_logical2=logical(loop_index); 
cal_data2=cal_data1 ; 
cal_data2(loopjndex_logical2)=0; 

%  find  largest  quiescent  period 
cal_data3=cal_data2; 

largest_quiescentJndex=ones(length(cal_data2),1); 

X=diff(LOCS); 

X=[LOCS(1)-1 ;  X;  length(cal_data2)-LOCS(length(LOCS))];  %  need  to  add  beginning  length 
l=find(X==max(X)); 


if  l==1  &&  length(time(round(length(time)*.05):LOCS(l)))>=60*60*Fs 
largest_quiescent_jndex(round(length(time)*.05):LOCS(l))=0; 
elseif  l==1  &&  length(time(round(length(time)*.05):LOCS(l)))<60*60*Fs 
l=find(X==max(X(2:length(X)))); 
end 

if  l==length(X)  &&  length(time(LOCS(l-1):length(cal_data2)-round(.05*length(time))))>=60*60*Fs 
largest_quiescentJndex(LOCS(l-1):length(cal_data2)-round(.05*length(time)))=0; 
elseif  l==length(X)  &&  length(time(LOCS(l-1):length(cal_data2)- 
round(.05*length(time))))<60*60*Fs 
l=find(X==max(X(2:length(X)-1))); 
end 

if  l>1  &&  klength(X); 

largest_quiescent_index(LOCS(l-1):LOCS(l))=0; 

end 

%  Ensure  X  is  at  least  an  hour 
Y  =  X(I); 
if  Y<=60*60*Fs 

error='quiscent  period  not  an  hour  use  high  amp  vibration  analysis' 
elseif  Y>60*60*Fs 
display('hour  of  quiescent  found') 
end 

largest_quiescentJndex=logical(largest_quiescentJndex); 

cal_data3(largest_quiescent_index)=0; 

channel1_cal_data=rsac_channel_1 ;  %channel  1 
channel1_cal_data(largest_quiescentjndex)=0; 

channel2_cal_data=rsac_channel_2;  %channel2 
channel2_cal_data(largest_quiescent_index)=0; 
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time_quiescent=time(~largest_quiescent_jndex); 

%  plot  longest  quiescent  period 
figure 

plot(time_quiescent,cal_data_qplot(~largest_quiescentJndex)) 
hold  on 

plot([time_quiescent(1 )  time_quiescent(length(time_quiescent))],[max_amp  max_amp],'r-') 
hold  on 

plot([time_quiescent(1 )  time_quiescent(length(time_quiescent))],[-max_amp  -max_amp],'r-') 
hold  on 

plot(time_quiescent,rsac1_qplot(~largest_quiescent_index),'g') 
hold  on 

plot(time_quiescent,rsac2_qplot(~largest_quiescentJndex),'m') 
ylabel('Velocity  (m/s)','Fontsize',12) 
xlabel('Time  (s)','Fontsize’,12) 

title('Longest  Continuous  Quiescent  Period  (Based  on  Channel  Z)',’Fontsize',16) 
legend('Calibrated  amplitude,  channel  Z max  amplitude', '+/-  max  amplitude'/Calibrated 
amplitude,  channel  1 '/Calibrated  amplitude,  channel  2') 

%  Power  Spectral  Density  of  quiescent  period 


N=2A1 2; 
overlap=N/2; 

figure 

[Pxx,F]=pwelch(cal_data3(~largest_quiescent_jndex),hamming(N), overlap, N,Fs); 

[Pxxl ,  FI  ]=pwelch(channel1_cal_data(~largest_quiescent_index), hamming(N), overlap, N,Fs); 
[Pxx2,F2]=pwelch(channel2_cal_data(~largest_quiescent_index), hamming(N), overlap, N,Fs); 
plot(F,1 0*log1 0(Pxx),'r') 
hold  on 

plot(F1 ,1 0*log1 0(Pxx1  ),'b') 
hold  on 

plot(F2,1 0*log1 0(Pxx2),'m') 

ylabel('dB  re  m/s/$$\sqrt{Hz}$$’, ’Interpreter', ’Latex’, 'Fontsize', 14) 
xlabel('Frequency  (Hz)','Fontsize',1 2) 

title({’Power  Spectral  Density  Estimate', 'Longest  Continuous  Quiescent  Period’}, 'Fontsize', 16) 
grid  on 

legend('Channel  Z', 'Channel  1 '/Channel  2') 
figure 

plot(F,1 0*log1 0(Pxx)/r’) 
hold  on 

plot(F1 ,1 0*log1 0(Pxx1  ),'b') 
hold  on 

plot(F2,1  CPIogl  0(Pxx2),’m') 

ylabel('dB  re  m/s/$$\sqrt{Hz}$$'/lnterpreter','Latex', 'Fontsize', 14) 
xlabel('Frequency  (Hz)’/Fontsize',1 2) 
xlim([0  2]) 

title({'Power  Spectral  Density  Estimate'/Longest  Continuous  Quiescent  Period,  0  to  2 
Hz'}, 'Fontsize', 16) 
grid  on 

legend(’Channel  Z’/Channel  1 '/Channel  2') 


length_z=length(cal_data); 
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taperz=floor(.05*length_z); 


lengthl  =length(rsac_channel_1 ); 
taperl  =floor(.05*length_1 ); 

length_2=length(rsac_channel_2); 

taper2=floor(.05*length_2); 

cal_data=cal_data(taperz:length_z-taperz); 
rsac_channel_1  =rsac_channel_1  (taperl  :length_1  -taperl ); 
rsac_channel_2=rsac_channel_2(taper2:length_2-taper2); 

time=time(taperz:length_z-taperz); 

%  Find  percentage  of  time  above  maximum  amplitude 

index=abs(cal_data)>=max_amp; 

above_max_amp_time=time(index); 

percentage_time_above_max_amp=(length(above_max_amp_time)/length(time))*100 

%  Create  square  wave  where  amplitude  surpasses  maximum  amplitude  for  use  in 
%  histogram. 

cal_data_update=cal_data; 
outlier=abs(cal_data_update)  >  max_amp; 

cal_data_update(~outiier)=0;  %  Sets  entire  signal  below  max  amp  to  zero 
cal_data_update(outlier)=sign(cal_data_update(outlier));  %  This  sets  all 
%  magnitudes  above  max  amp  to  1  (builds  a  square  wave) 

W=pulsewidth(abs(cal_data_update),Fs);  %  stores  times  of  pulses  above  max  amp 

%  W  =  pulsewidth(X)  returns  a  vector,  W,  containing  the  time  differences 
%  between  the  mid-reference  level  instants  of  the  initial  and  final  transitions 
%  of  each  positive-polarity  pulse  in  the  bilevel  waveform,  X. 

[S,INITCROSS]=pulsesep(abs(cal_data_update),Fs); 

S=[INITCF!OSS(1);  S]; 

%  S  =  pulsesep(X)  returns  the  differences,  S,  between  the  mid-reference  level 
%instants  of  the  final  negative-going  transitions  of  every  positive-polarity 
%pulse  and  the  next  positive-going  transition.  X  is  a  bilevel  waveform. 

%  [SJNITCROSS]  =  pulsesep(...)  returns  the  mid-reference  level  instants, 

%  INITCROSS,  of  the  first  positive-polarity  transitions. 

figure 

histogram(W) 

set(gca,  Yscale','log')  %  this  gives  you  average  time  of  pulses  above  max  amp 
title({'Histogram  of  Length  of  Continuous  Time  Signal  Exceeding  Maximum', 'Channel 
Z'},'Fontsize',16) 
xlabel('Time  (s)',’Fontsize',12) 

figure 

histogram(W) 
ylim([0  50]) 
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xlabel('Time  (s)','Fontsize',12) 

title({’Histogram  of  Length  of  Continuous  Time  Signal  Exceeding  Maximum', 'Channel 
Z'},'Fontsize',16) 

figure 

histogram(S) 
set(gca,  Yscale','log') 
xlabel('Time  (s)','Fontsize’,12) 

title({’Flistogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude’, ’Channel 
Z'},'Fontsize’,16) 

figure 

histogram(S) 
ylim([0  50]) 

xlabel('Time  (s)','Fontsize’,12) 

titlett'Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude’, ’Channel 
Z'},'Fontsize’,16) 

figure 

histogram(cal_data,400) 
set(gca,  Yscale',’log') 
xlabel('Velocity  Amplitude  (m/s)') 
ylabel('Number  of  Occurrences') 
title(’Histogram  of  Channel  Z  Amplitudes') 

figure 

histogram(cal_data,400) 
xlabel('Velocity  Amplitude  (m/s)') 
ylabel('Number  of  Occurrences  maxed  at  20') 
ylim([0  20]) 

title('Histogram  of  Channel  Z  Amplitudes') 
std_z=std(cal_data) 

B.  HIGH  AMPLITUDE  VIBRATION  ANALYSIS  CODE 

function 

high_amp_vibration_analysis(rsac_channel_Z,max_amp,Fs,rsac_channel_1  ,rsac_channel_2) 
%  OBS,  ocean  bottom  seismometer,  maximum  amplitude,  BHZ,  sampling  frequency 

% 

%  Purpose:  Run  this  function  if  the  sample  day  too  noisy  for 
%  vibration_analysis.m  (not  at  least  1  hour  in  length  where  the  amplitude 
%  does  not  exceed  maximum  amplitude).  Takes  three  calibrated  channels  of 
%  OBS  data,  OBS  sampling  frequency,  and  maximum  desired  amplitude,  and 
%  finds  the  three  hours  with  lowest  mean  amplitude.  Function  outputs  an 
%  amplitude  (m/s)  vs  time  plot  (for  all  three  channels),  the  mean  and 
%  standard  deviation  of  entire  sample  amplitude,  power  spectral  densities 
%  of  three-hour  with  lowest  mean,  a  percentage  of  time  spent  above  maximum 
%  desired  amplitude,  a  histogram  of  amplitudes,  a  histogram  of  lengths  of 
%  signals  greater  than  max  desired  amplitude,  and  a  histogram  of  time 
%  periods  between  signals  that  exceed  maximum  amplitude. 

% 

%  Sample  needs  to  be  at  least  1 5  hours  long  for  this  function  to  work 
%  properly. 
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% 

%  Output  Variables: 

%  std_z  =  standard  deviation  of  channel  z  amplitude 
%  percentage_time_above_max_amp  =  percentage  of  time  where  amplitude 
%  is  greater  than  max  amplitude 

% 

%  Input  Variables: 

%  rsac_channel_Z  =  vertical  .sac  file  after  it  has  been  read  through 
%  rsac.m 

%  max_amp  =  desired  maximum  amplitude  (m/s) 

%  Fs  =  sampling  frequency  of  OBS  (Hz) 

%  rsac_channel_1  =  channel  1  .sac  file  after  it  has  been  read  through 
%  rsac.m 

%  rsac_channel_2  =  channel  2  .sac  file  after  it  has  been  read  through 
%  rsac.m 

% 

%  Local  Variables: 

%  cal_data  =  calibrated  vertical  channel  amplitude  (m/s) 

%  time  =  time  (seconds) 

%  length_z  =  length  of  cal_data  (vector  length),  channel  Z 
%  taperz  =  amount  to  subtract  from  beginning/end  of  sample  to  get  rid 
%  of  taper,  channel  Z 

%  length_1  =  vector  length  of  channel  1  amplitude  vector 
%  taperl  =  amount  to  subtract  from  beginning/end  of  sample  to  get  rid 
%  of  taper,  channel  1 

%  length_2  =  vector  length  of  channel  2  amplitude  vector 
%  taper2  =  amount  to  subtract  from  beginning/end  of  sample  to  get  rid 
%  of  taper,  channel  2 

%  rsac_channel_1  =  amplitude  vector  channel  1 

%  rsac_channel_2  =  amplitdue  vector  channel  2 

%  second_3_hours_cal_data  to  twentyone_3_hours_cal_data  =  three-hour 
%  blocks  of  time  used  to  find  three-hour  time  period  with  lowest  mean 
%  ampltidue 

%  MEAN  =  vector  containing  mean  values  of  three-hour  blocks  of  time 

%  Y  =  used  in  first  if/elseif  statement.  If  size  of  sample  much  less 

%  than  24  hours,  Y  variable  helps  assign  correct  area  to  apply  taper 
%  X  =  used  in  if/elseif  statement  to  choose  correct  three-hour  block 
%  and  assign  variables  to  PSD 

%  N  =  FFT  length 

%  overlap  =  pwelch  overlap  value 

%  Pxx  =  channel  Z  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  F,  FI ,  F2  =  Frequency  (Hz) 

%  Pxxl  =  channel  1  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  Pxx2  =  channel  2  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  index  =  index  locating  where  cal_data  is  greater  than  the  maximum 

%  amplitude 

%  above_max_amp_time  =  time  where  cal_data  is  greater  than  maximum 
%  amplitude  (s) 

%  outlier  =  logical  index  location  where  Z-amplitude  is  greater  than 
%  maximum  amplitude 

%  cal_data_update  =  same  as  cal_data,  but  only  contains  amplitudes 
%  greater  than  maximum  amplitude,  (m/s) 
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%  W  =  vector  with  lengths  of  signals  that  exceed  max  amplitude 
%  S  =  vector  that  contains  different  lengths  of  time  between  signals 
%  that  exceed  max  amplitude 

% 

%  Functions  Called:  None. 

% 

%  Filename:  high_amp_vibration_analysis.m 

%  Written  by:  Jeremy  R.  Hankins 

cal_data=rsac_channel_Z(:,2); 

time=rsac_channel_Z(:,1); 

length_z=length(cal_data); 

taperz=floor(.05*length_z); 

lengthl  =length(rsac_channel_1  (:,2)); 
taperl  =floor(.05*length_1 ); 

length_2=length(rsac_channel_2(:,2)); 

taper2=floor(.05*length_2); 

rsac_channel_1=rsac_channel_1(:,2); 

rsac_channel_2=rsac_channel_2(:,2); 

figure 

plot(time,cal_data) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],’r-') 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
xlabel('Time  (s)','Fontsize’,12) 

title(’Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Vertical  Velocity  (Channel 
Z)7Fontsize',16) 

legend('Calibrated  amplitude', ’+/-  max  amplitude') 
figure 

plot(time,rsac_channel_1 ) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],’r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],’r-') 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
xlabel('Time  (s)','Fontsize',12) 

title('Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity 

(Channel  1)','Fontsize',16) 

legend('Calibrated  amplitude', '+/-  max  amplitude’) 

figure 

plot(time,rsac_channel_2) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],’r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],'r~') 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
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xlabel('Time  (s)','Fontsize',12) 

title('Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity 

(Channel  2)','Fontsize’,16) 

legend('Calibrated  amplitude','+/-  max  amplitude’) 

%  Divide  sample  into  3  hour  time  blocks,  exluding  taper  at  beginning  and 
%  end  of  sample.  3  hour  block  that  includes  taper  could  be  less/more 
%  than  three  hours.  Amount  less/more  depends  on  length  of  sample. 

%  if/elseif  statement  covers  files  as  short  as  15  hours 

second_3_hours_cal_data=cal_data(taperz:4*60*60*Fs); 
third_3_hours_cal_data=cal_data(2*60*60*Fs:5*60*60*Fs); 
fourth_3_hours_cal_data=cal_data(3*60*60*Fs:6*60*60*Fs); 
fifth_3_hours_cal_data=cal_data(4*60*60*Fs:7*60*60*Fs); 
sixth_3_hours_cal_data=cal_data(5*60*60*Fs:8*60*60*Fs); 
seventh_3_hours_cal_data=cal_data(6*60*60*Fs:9*60*60*Fs); 
eighth_3_hours_cal_data=cal_data(7*60*60*Fs:10*60*60*Fs); 
nine_3_hours_cal_data=cal_data(8*60*60*Fs:1 1*60*60*Fs); 
ten_3_hours_cal_data=cal_data(9*60*60*Fs:12*60*60*Fs); 
eleven_3_hours_cal_data=cal_data(1 0*60*60*Fs:1 3*60*60*Fs); 
twelve_3_hours_cal_data=cal_data(1 1  *60*60*Fs:1 4*60*60*Fs); 
if  length(cal_data)/(Fs*60*60)>=23; 

Y=20; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(1 5*60*60*Fs:1 8*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(1 6*60*60*Fs:1 9*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:20*60*60*Fs); 
nineteen_3_hours_cal_data=cal_data(1 8*60*60*Fs:21  *60*60*Fs); 
twenty_3_hours_cal_data=cal_data(19*60*60*Fs:22*60*60*Fs); 
twentyone_3_hours_cal_data=cal_data(20*60*60*Fs:length_z-taperz); 
elseif  length(cal_data)/(Fs*60*60)>=22  &&  length(cal_data)/(Fs*60*60)  <  23; 

Y=19; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(1 5*60*60*Fs:1 8*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(1 6*60*60*Fs:1 9*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:20*60*60*Fs); 
nineteen_3_hours_cal_data=cal_data(1 8*60*60*Fs:21  *60*60*Fs); 
twenty_3_hours_cal_data=cal_data(19*60*60*Fs:length_z-taperz); 
twentyone_3_hours_cal_data=1  El  0; 

elseif  length(cal_data)/(Fs*60*60)>=21  &&  length(cal_data)/(Fs*60*60)  <  22; 

Y=18; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(1 5*60*60*Fs:1 8*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(1 6*60*60*Fs:1 9*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:20*60*60*Fs); 
nineteen_3_hours_cal_data=cal_data(18*60*60*Fs:length_z-taperz); 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 
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elseif  length(cal_data)/(Fs*60*60)>=20  &&  length(cal_data)/(Fs*60*60)  <  21 ; 
Y=17; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(1 5*60*60*Fs:1 8*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(1 6*60*60*Fs:1 9*60*60*Fs); 
eighteen_3_hours_cal_data=cal_data(17*60*60*Fs:length_z-taperz); 
nineteen_3_hours_cal_data=1  El  0; 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 

elseif  length(cal_data)/(Fs*60*60)>=19  &&  length(cal_data)/(Fs*60*60)  <  20; 
Y=16; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(1 5*60*60*Fs:1 8*60*60*Fs); 
seventeen_3_hours_cal_data=cal_data(16*60*60*Fs:length_z-taperz); 
eighteen_3_hours_cal_data=1  El  0; 
nineteen_3_hours_cal_data=1  El  0; 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 

elseif  length(cal_data)/(Fs*60*60)>=18  &&  length(cal_data)/(Fs*60*60)  <  19; 

Y=15; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:17*60*60*Fs); 
sixteen_3_hours_cal_data=cal_data(15*60*60*Fs:length_z-taperz); 
seventeen_3_hours_cal_data=1  El  0; 
eighteen_3_hours_cal_data=1  El  0; 
nineteen_3_hours_cal_data=1  El  0; 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 

elseif  length(cal_data)/(Fs*60*60)>=17  &&  length(cal_data)/(Fs*60*60)  <  18; 
Y=14; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(1 3*60*60*Fs:1 6*60*60*Fs); 
fifteen_3_hours_cal_data=cal_data(14*60*60*Fs:length_z-taperz); 
sixteen_3_hours_cal_data=1  El  0; 
seventeen_3_hours_cal_data=1  El  0; 
eighteen_3_hours_cal_data=1  El  0; 
nineteen_3_hours_cal_data=1  El  0; 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 

elseif  length(cal_data)/(Fs*60*60)>=16  &&  length(cal_data)/(Fs*60*60)  <17; 

Y=13; 

thirteen_3_hours_cal_data=cal_data(1 2*60*60*Fs:1 5*60*60*Fs); 
fourteen_3_hours_cal_data=cal_data(13*60*60*Fs:length_z-taperz); 
fifteen_3_hours_cal_data=1  El  0; 
sixteen_3_hours_cal_data=1  El  0; 
seventeen_3_hours_cal_data=1  El  0; 
eighteen_3_hours_cal_data=1  El  0; 
nineteen_3_hours_cal_data=1  El  0; 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 
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elseif  length(cal_data)/(Fs*60*60)>=15  &&  length(cal_data)/(Fs*60*60)  <  16; 

Y=1 2; 

thirteen_3_hours_cal_data=cal_data(12*60*60*Fs:length_z-taperz); 
fourteen_3_hours_cal_data=1  El  0; 
fifteen_3_hours_cal_data=1  El  0; 
sixteen_3_hours_cal_data=1  El  0; 
seventeen_3_hours_cal_data=1  El  0; 
eighteen_3_hours_cal_data=1  El  0; 
nineteen_3_hours_cal_data=1  El  0; 
twenty_3_hours_cal_data=1  El  0; 
twentyone_3_hours_cal_data=1  El  0; 
end 

MEAN=[mean(abs(second_3_hours_cal_data)),... 

mean(abs(third_3_hours_cal_data)),mean(abs(fourth_3_hours_cal_data)),... 

mean(abs(fifth_3_hours_cal_data)),mean(abs(sixth_3_hours_cal_data)),... 

mean(abs(seventh_3_hours_cal_data)),mean(abs(eighth_3_hours_cal_data)),... 

mean(abs(nine_3_hours_cal_data)),mean(abs(ten_3_hours_cal_data)),... 

mean(abs(eleven_3_hours_cal_data)),mean(abs(twelve_3_hours_cal_data)),... 

mean(abs(thirteen_3_hours_cal_data)),mean(abs(fourteen_3_hours_cal_data)),... 

mean(abs(fifteen_3_hours_cal_data)),mean(abs(sixteen_3_hours_cal_data)),... 

mean(abs(seventeen_3_hours_cal_data)),mean(abs(eighteen_3_hours_cal_data)),. 

mean(abs(nineteen_3_hours_cal_data)),mean(abs(twenty_3_hours_cal_data)),... 

mean(abs(twentyone_3_hours_cal_data))]; 

X=find(MEAN==min(MEAN)) 

%  if/elseif  statments  to  plot  correct  psd 


if  X==1 

psd_z=second_3_hours_cal_data; 
psd_1  =rsac_channel_1  (taperl  :4*60*60*Fs); 
psd_2=rsac_channel_2(taper2:4*60*60*Fs); 
elseif  X==2 

psd_z=third_3_hours_cal_data; 
psd_1  =rsac_channel_1  (2*60*60*Fs:5*60*60*Fs); 
psd_2=rsac_channel_2(2*60*60*Fs:5*60*60*Fs); 
elseif  X==3 

psd_z=fourth_3_hours_cal_data; 
psd_1  =rsac_channel_1  (3*60*60*Fs:6*60*60*Fs); 
psd_2=rsac_channel_2(3*60*60*Fs:6*60*60*Fs); 
elseif  X==4 

psd_z=fifth_3_hours_cal_data; 
psd_1  =rsac_channel_1  (4*60*60*Fs:7*60*60*Fs); 
psd_2=rsac_channel_2(4*60*60*Fs:7*60*60*Fs); 
elseif  X==5 

psd_z=sixth_3_hours_cal_data; 
psd_1  =rsac_channel_1  (5*60*60*Fs:8*60*60*Fs); 
psd_2=rsac_channel_2(5*60*60*Fs:8*60*60*Fs); 
elseif  X==6 

psd_z=seventh_3_hours_cal_data; 
psd_1  =rsac_channel_1  (6*60*60*Fs:9*60*60*Fs); 
psd_2=rsac_channel_2(6*60*60*Fs:9*60*60*Fs); 
elseif  X==7 

psd_z=eighth_3_hours_cal_data; 
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pscM  =rsac_channel_1  (7*60*60*Fs:1 0*60*60*Fs); 
psd_2=rsac_channel_2(7*60*60*Fs:10*60*60*Fs); 
elseif  X==8 

psd_z= n  i  n  e_3_h  o  u  rs_cal_d  ata ; 
psd_1  =rsac_channel_1  (8*60*60*Fs:1 1  *60*60*Fs); 
psd_2=rsac_channel_2(8*60*60*Fs:1 1  *60*60*  Fs); 
elseif  X==9 

psd_z=te  n_3_h  o  u  rs_cal_d  ata ; 
psd_1  =rsac_channel_1  (9*60*60*Fs:1 2*60*60*Fs); 
psd_2=rsac_channel_2(9*60*60*Fs:12*60*60*Fs); 
elseif  X==1 0 

psd_z=eleven_3_hours_cal_data; 
psd_1  =rsac_channel_1  (1 0*60*60*Fs:1 3*60*60*Fs); 
psd_2=rsac_channel_2(1 0*60*60*Fs:1 3*60*60*Fs); 
elseif  X==1 1 

psd_z=twelve_3_hours_cal_data; 
psd_1=rsac_channel_1  (1 1*60*60*Fs:14*60*60*Fs); 
psd_2=rsac_channel_2(1 1*60*60*Fs:14*60*60*Fs); 
elseif  X==1 2 

psd_z=thirteen_3_hours_cal_data; 
if  Y==1 2 

psd_1  =rsac_channel_1  (1 2*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(12*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 2*60*60*Fs:1 5*60*60*Fs); 
psd_2=rsac_channel_2(1 2*60*60*Fs:1 5*60*60*Fs); 
end 

elseif  X==13 

psd_z=fourteen_3_hours_cal_data; 
if  Y==1 3 

psd_1  =rsac_channel_1  (1 3*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(13*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 3*60*60*Fs:1 6*60*60*Fs); 
psd_2=rsac_channel_2(1 3*60*60*Fs:1 6*60*60*Fs); 
end 

elseif  X==14 

psd_z=fifteen_3_hours_cal_data; 
if  Y==14 

psd_1  =rsac_channel_1  (1 4*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(14*60*60*Fs:length_2-taper2); 
else 

psd_1=rsac_channel_1(14*60*60*Fs:1 7*60*60*Fs); 
psd_2=rsac_channel_2(14*60*60*Fs:1 7*60*60*Fs); 
end 

elseif  X==15 

psd_z=sixteen_3_hours_cal_data; 
if  Y==1 5 

psd_1  =rsac_channel_1  (1 5*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(15*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 5*60*60*Fs:1 8*60*60*Fs); 
psd_2=rsac_channel_2(1 5*60*60*Fs:1 8*60*60*Fs); 
end 

elseif  X==1 6 
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psd_z=seventeen_3_hours_cal_data; 
if  Y==1 6 

psd_1  =rsac_channel_1  (1 6*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(16*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 6*60*60*Fs:1 9*60*60*Fs); 
psd_2=rsac_channel_2(1 6*60*60*Fs:1 9*60*60*Fs); 
end 

elseif  X==1 7 

psd_z=eighteen_3_hours_cal_data; 
if  Y==17 

psd_1  =rsac_channel_1  (1 7*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(17*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 7*60*60*Fs:20*60*60*Fs); 
psd_2=rsac_channel_2(17*60*60*Fs:20*60*60*Fs); 
end 

elseif  X==18 

psd_z=nineteen_3_hours_cal_data; 
if  Y==1 8 

psd_1  =rsac_channel_1  (1 8*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(18*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 8*60*60*Fs:21  *60*60*Fs); 
psd_2=rsac_channel_2(1 8*60*60*Fs:21  *60*60*Fs); 
end 

elseif  X==19 

psd_z=twenty_3_hours_cal_data; 
if  Y==19; 

psd_1  =rsac_channel_1  (1 9*60*60*Fs:length_1  -taperl ); 
psd_2=rsac_channel_2(19*60*60*Fs:length_2-taper2); 
else 

psd_1  =rsac_channel_1  (1 9*60*60*Fs:22*60*60*Fs); 
psd_2=rsac_channel_2(19*60*60*Fs:22*60*60*Fs); 
end 

elseif  X==20 

psd_z=twentyone_3_hours_cal_data; 

psd_1  =rsac_channeM  (20*60*60*Fs:length_1  -taperl ); 

psd_2=rsac_channel_2(20*60*60*Fs:length_2-taper2); 

end 


N=2A1 2; 
overlap=N/2; 

figure 

[Pxx,F]=pwelch(psd_z,hamming(N), overlap, N,Fs); 

[Pxxl  ,F1]=pwelch(psd_1  ,hamming(N),overlap,N,Fs); 
[Pxx2,F2]=pwelch(psd_2,hamming(N),overlap,N,Fs); 
plot(F,1 0*log1 0(Pxx),'r') 
hold  on 

plot(F1 ,1 0*log1 0(Pxx1  ),'b') 
hold  on 

plot(F2,1 0*log1 0(Pxx2),'m') 

ylabel('dB  re  m/s/$$\sqrt{Hz}$$','lnterpreterYLatexYFontsize',14) 
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xlabel('Frequency  (Hz)’,'Fontsize',1 2) 
xlim([0  2]) 

title({'Power  Spectral  Density  Estimate', 'Three  Flours  with  Lowest  Mean  Amplitude,  0  to  2 

Hz’},'Fontsize',16) 

grid  on 

legend('Channel  Z’, ’Channel  1  '/Channel  2') 
figure 

plot(F,10*log10(Pxx),'r') 
hold  on 

plot(F1,10*log10(Pxx1),'b') 
hold  on 

plot(F2,10*log10(Pxx2),’m') 

ylabel(’dB  re  m/s/$$\sqrt{Hz}$$’, ’Interpreter’, ’Latex’, 'Fontsize', 14) 
xlabel('Frequency  (Hz)', ’Fontsize', 1 2) 

title({'Power  Spectral  Density  Estimate’, ’Three  Hours  with  Lowest  Mean  Amplitude’},’Fontsize',16) 
grid  on 

legend(’Channel  Z’, ’Channel  1’, ’Channel  2’) 

%  get  rid  of  taper  for  plot  and  histograms 
cal_data=cal_data(taperz:length_z-taperz); 
rsac_channel_1  =rsac_channel_1  (taperl  :length_1  -taperl ); 
rsac_channel_2=rsac_channel_2(taper2:length_2-taper2); 
time=time(taperz:length_z-taperz); 

%  psd  for  entire  day  of  data 
figure 

[Pxx,F]=pwelch(cal_data,hamming(N),  overlap,  N,Fs); 

[Pxxl , FI  ]=pwelch(rsac_channel_1,hamming(N), overlap, N,Fs); 
[Pxx2,F2]=pwelch(rsac_channel_2,hamming(N), overlap, N,Fs); 
plot(F,10*log10(Pxx),’r') 
hold  on 

plot(F1 ,1 0*log1 0(Pxx1  ),’b') 
hold  on 

plot(F2,10*log10(Pxx2),’m') 

ylabel(’dB  re  m/s/$$\sqrt{Hz}$$’, ’Interpreter’, ’Latex’, 'Fontsize', 14) 
xlabel(’Frequency  (Hz)’, 'Fontsize', 1 2) 

title({'Power  Spectral  Density  Estimate’, ’Entire  Day  of  Data’}, 'Fontsize', 16) 
grid  on 

legend(’Channel  Z’, ’Channel  1’, ’Channel  2’) 

index=abs(cal_data)>=max_amp; 

above_max_amp_time=time(index); 

percentage_time_above_max_amp=(length(above_max_amp_time)/length(time))*100 

cal_data_update=cal_data; 
outlier=abs(cal_data_update)  >  max_amp; 

cal_data_update(~outiier)=0;  %  Sets  entire  signal  below  max  amp  to  zero 

cal_data_update(outlier)=sign(cal_data_update(outlier));  %  This  sets  all  magnitudes  above  max 
amp  to  1  (builds  a  square  wave) 

W=pulsewidth(abs(cal_data_update),Fs);  %  stores  times  of  pulses  above  max  amp 

[S,INITCROSS]=pulsesep(abs(cal_data_update),Fs); 

S=[INITCROSS(1);  S]; 
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figure 

histogram(W) 

set(gca,  Yscale','log')  %  this  gives  you  average  time  of  pulses  above  max  amp 
title({'Histogram  of  Length  of  Continuous  Time  Signal  Exceeding  Maximum', 'Channel 
Z'},'Fontsize',16) 
xlabel('Time  (s)’,'Fontsize’,12) 

figure 

histogram(W) 
ylim([0  50]) 

xlabel('Time  (s)','Fontsize’,12) 

title({’Histogram  of  Length  of  Continuous  Time  Signal  Exceeding  Maximum’, 'Channel 
Z'},’Fontsize’,16) 

figure 

histogram(S) 
set(gca,  YscaleVlog') 
xlabel('Time  (s)','Fontsize',12) 

title({'Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude', 'Channel 
Z'},'Fontsize',16) 

figure 

histogram(S) 
ylim([0  50]) 

xlabel('Time  (s)','Fontsize',12) 

title({’Histogram  of  Length  of  Time  Between  Signals  that  Exceed  Max  Amplitude', 'Channel 
Z'},'Fontsize',16) 

figure 

histogram(cal_data,400) 
set(gca,  YscaleVlog') 
xlabel('Velocity  Amplitudes  (m/s)') 
ylabel('Number  of  Occurrences') 
title('Histogram  of  Channel  Z  Amplitudes') 

figure 

histogram(cal_data,400) 
xlabel('Velocity  Amplitudes  (m/s)') 
ylabel('Number  of  Occurrences  maxed  at  20') 
ylim([0  20]) 

title(’Histogram  of  Channel  Z  Amplitudes') 
std_z=std(cal_data) 
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c. 


LOW  AMPLITUDE  VIBRATION  ANALYSIS  CODE 


function 

low_amp_vibration_analysis(rsac_channel_Z,max_amp,Fs,rsac_channel_1  ,rsac_channel_2) 
%  OBS,  ocean  bottom  seismometer,  maximum  amplitude,  BHZ,  sampling  frequency 

% 

%  Purpose:  Run  this  function  if  the  sample  day  does  not  exceed  max 
%  amplitude  with  vibration_analysis.m  (does  not  cross  max  amplitude  at 
%  all).  Takes  three  calibrated  channels  of  OBS  data,  OBS  sampling 
%  frequency,  and  maximum  desired  amplitude.  Function  outputs  an 
%  amplitude  (m/s)  vs  time  plot  (for  all  three  channels),  the  mean  and 
%  standard  deviation  of  entire  sample  amplitude,  power  spectral  densities 
%  of  entire  day,  and  a  histogram  of  amplitudes. 

% 

%  Sample  needs  to  be  at  least  1 5  hours  long  for  this  function  to  work 
%  properly. 

% 

%  Output  Variables: 

%  avg_z  =  mean  of  absolute  value  of  channel  z  amplitude 
%  std_z  =  standard  deviation  of  absolute  value  of  channel  z  amplitude 
%  percentage_time_above_max_amp  =  percentage  of  time  above  max  amp 
%  (should  be  zero) 

% 

%  Input  Variables: 

%  rsac_channel_Z  =  vertical  .sac  file  after  it  has  been  read  through 
%  rsac.m 

%  max_amp  =  desired  maximum  amplitude  (m/s) 

%  Fs  =  sampling  frequency  of  OBS  (Hz) 

%  rsac_channeM  =  channel  1  .sac  file  after  it  has  been  read  through 
%  rsac.m 

%  rsac_channel_2  =  channel  2  .sac  file  after  it  has  been  read  through 
%  rsac.m 
% 

%  Local  Variables: 

%  cal_data  =  calibrated  vertical  channel  amplitude  (m/s) 

%  time  =  time  (seconds) 

%  length_z  =  length  of  cal_data  (vector  length),  channel  Z 
%  taperz  =  amount  to  subtract  from  beginning/end  of  sample  to  get  rid 
%  of  taper,  channel  Z 

%  length_1  =  vector  length  of  channel  1  amplitude  vector 
%  taperl  =  amount  to  subtract  from  beginning/end  of  sample  to  get  rid 
%  of  taper,  channel  1 

%  length_2  =  vector  length  of  channel  2  amplitude  vector 
%  taper2  =  amount  to  subtract  from  beginning/end  of  sample  to  get  rid 
%  of  taper,  channel  2 

%  rsac_channel_1  =  amplitude  vector  channel  1 

%  rsac_channel_2  =  amplitdue  vector  channel  2 

%  N  =  FFT  length 
%  overlap  =  pwelch  overlap  value 

%  Pxx  =  channel  Z  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  F,  FI ,  F2  =  Frequency  (Hz) 

%  Pxxl  =  channel  1  distribution  of  power  per  unit  of  freq,  but  in  this 
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%  application  units  of  [m/s/sqrt(Hz)] 

%  Pxx2  =  channel  2  distribution  of  power  per  unit  of  freq,  but  in  this 
%  application  units  of  [m/s/sqrt(Hz)] 

%  index  =  index  locating  where  cal_data  is  greater  than  the  maximum 
%  amplitude 

%  above_max_amp_time  =  time  where  cal_data  is  greater  than  maximum 
%  amplitude  (s) 

% 

%  Functions  Called:  None. 

% 

%  Filename:  low_amp_vibration_analysis.m 

%  Written  by:  Jeremy  R.  Hankins 

cal_data=rsac_channel_Z(:,2); 

time=rsac_channel_Z(:,1); 

length_z=length(cal_data); 

taperz=floor(.05*length_z); 

length  1  =length(rsac_channel_1  (:,2)); 
taperl  =floor(.05*length_1 ); 

length_2=length(rsac_channel_2(:,2)); 

taper2=floor(.05*length_2); 

rsac_channel_1=rsac_channel_1(:,2); 

rsac_channel_2=rsac_channel_2(:,2); 

figure 

plot(time,cal_data) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],'r-') 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
xlabel('Time  (s)','Fontsize’,12) 

title(’Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Vertical  Velocity  (Channel 
Z)YFontsize\16) 

legend('Calibrated  amplitude', ’+/-  max  amplitude') 
figure 

plot(time,rsac_channel_1 ) 
hold  on 

plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp],'r-’) 
ylabel('Velocity  (m/s)','Fontsize',1 2) 
xlabel('Time  (s)','Fontsize',12) 

title('Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity 

(Channel  1)','Fontsize',16) 

legend('Calibrated  amplitude', '+/-  max  amplitude’) 

figure 

plot(time,rsac_channel_2) 
hold  on 
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plot([time(1)  time(length(time))],[max_amp  max_amp],'r-') 
hold  on 

plot([time(1)  time(length(time))],[-max_amp  -max_amp]/r-’) 
ylabel('Velocity  (m/s)',’Fontsize',1 2) 
xlabel('Time  (s)','Fontsize',12) 

title('Entire  Day  of  Data  with  Desired  Maximum  Amplitude  Displayed,  Horizontal  Velocity 

(Channel  2)', 'Fontsize', 16) 

legend('Calibrated  amplitude’, ’+/-  max  amplitude’) 

%  get  rid  of  tapers 

cal_data=cal_data(taperz:length_z-taperz); 
rsac_channel_1  =rsac_channel_1  (taperl  :length_1  -taperl ); 
rsac_channel_2=rsac_channel_2(taper2:length_2-taper2); 
time=time(taperz:length_z-taperz); 


N=2A1 2; 
overlap=N/2; 

figure 

[Pxx,F]=pwelch(cal_data,hamming(N), overlap,  N,Fs); 

[Pxxl ,  FI  ]=pwelch(rsac_channel_1,hamming(N),  overlap,  N,Fs); 
[Pxx2,F2]=pwelch(rsac_channel_2,hamming(N), overlap, N,Fs); 
plot(F,1  CPIogl  0(Pxx),’r') 
hold  on 

plot(F1 ,1 0*log1 0(Pxx1  ),'b') 
hold  on 

plot(F2,1  CPIogl  0(Pxx2),’m') 

ylabel(’dB  re  m/s/$$\sqrt{Hz}$$’, ’Interpreter’, ’Latex’, 'Fontsize', 14) 
xlabel(’Frequency  (Hz)’,'Fontsize',1 2) 
title({'Power  Spectral  Density  Estimate'}, ’Fontsize’, 16) 
grid  on 

legend(’Channel  Z’, ’Channel  1 ’/Channel  2') 
figure 

plot(F,1  CPIogl  0(Pxx),’r') 
hold  on 

plot(F1 ,1 0*log1 0(Pxx1  ),'b') 
hold  on 

plot(F2,1 0*log1 0(Pxx2),’m') 

ylabel('dB  re  m/s/$$\sqrt{Hz}$$’, ’Interpreter', 'Latex', 'Fontsize', 14) 
xlabel(’Frequency  (Hz)','Fontsize’,1 2) 
xlim([0  2]) 

title({’Power  Spectral  Density  Estimate', '0  to  2  Hz'}, 'Fontsize', 16) 
grid  on 

legend(’Channel  Z’, ’Channel  1 ’/Channel  2’) 

index=abs(cal_data)>=max_amp; 

above_max_amp_time=time(index); 

percentage_above_max_amp_time=(length(above_max_amp_time)/length(time))*100  %  will  be 
zero  with  this  function 

figure 

histogram(cal_data,400) 

set(gca/Yscale’/log’) 
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APPENDIX  B.  DATA  TABLES 


Table  3.  YM  39  and  YM  38  Standard  Deviation  and  Percentage  of 
Time  Spent  Above  Maximum  Amplitude  Threshold. 


YM  39 

YM38 

Date 

Standard 
deviation  (pm/s) 

RMS  acceleration 
(pm/s2) 

Percentage  above 
max  amplitude 

Date 

Standard 
deviation  (pm/s) 

RMS  acceleration 
(pm/s2) 

Percentage  above 
maxamplitude 

13-Jun-08 

191 

155 

17 

13-Jun-08 

0.23 

2.7 

0.00 

14-Jun-08 

143 

65.8 

10 

14-Jun-08 

0.40 

3.6 

0.00 

l-Jul-09 

112 

57.0 

36 

l-Jul-09 

0.45 

3.6 

0.00 

15-Jul-08 

604 

919 

83 

15-Jul-08 

1.9 

9.5 

0.03 

l-Aug-08 

395 

355 

86 

l-Aug-08 

7.3 

8.1 

0.16 

15-Aug-08 

324 

403 

83 

15-Aug-08 

1.2 

3.7 

0.00 

l-Sep-08 

662 

574 

91 

l-Sep-08 

0.57 

3.8 

0.00 

15-Sep-08 

350 

392 

86 

15-Sep-08 

0.70 

6.1 

0.00 

l-Oct-08 

365 

364 

83 

l-Oct-08 

0.56 

3.9 

0.00 

15-Oct-08 

375 

303 

77 

15-Oct-08 

0.44 

4.8 

0.00 

l-Nov-08 

399 

358 

80 

l-Nov-08 

1.1 

3.8 

0.01 

15-NOV-08 

219 

279 

74 

15-NOV-08 

1.2 

12.7 

0.00 

13-Dec-08 

200 

489 

42 

13-Dec-08 

1.3 

4.1 

0.01 

14-Dec-08 

277 

303 

33 

14-Dec-08 

1.9 

4.3 

0.04 

l-Jan-09 

727 

185 

84 

l-Jan-09 

0.55 

4.6 

0.00 

15-Jan-09 

366 

323 

68 

15-Jan-09 

1.3 

3.9 

0.02 

l-Feb-09 

46.4 

128 

10 

l-Feb-09 

0.42 

3.4 

0.00 

15-Feb-09 

150 

33.0 

30 

15-Feb-09 

0.43 

3.9 

0.00 

l-Mar-09 

309 

48.4 

57 

l-Mar-09 

1.3 

5.1 

0.00 

15-Mar-09 

119 

74.6 

54 

15-Mar-09 

0.57 

4.7 

0.00 

l-Apr-09 

96.9 

19.0 

34 

l-Apr-09 

2.9 

13.3 

0.03 

15-Apr-09 

88.0 

71.9 

8 

15-Apr-09 

1.3 

8.8 

0.01 

1-M  ay-09 

63.9 

86.3 

25 

l-May-09 

0.56 

4.7 

0.00 

15-May-09 

16.1 

28.9 

1 

15-May-09 

0.31 

3.4 

0.00 
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Table  4.  7D  J65A  and  7D  J73A  Standard  Deviation  and  Percentage  of 
Time  Spent  Above  Maximum  Amplitude  Threshold. 


7DJ65A 

7DJ73A 

Date 

Standard 
deviation  (pm/s) 

RMS  acceleration 
(pm/s2) 

Percentage  above 
max  amplitude 

Date 

Standard 
deviation  (pm/s) 

RMS  acceleration 
(pm/s2) 

Percentage  above 
max  amplitude 

1-Nov-ll 

30.5 

51.7 

9.8 

1-Nov-ll 

28.7 

33.0 

8 

15-Nov-ll 

26.8 

38.5 

6.6 

15-Nov-ll 

24.4 

25.2 

5 

1-Dec-ll 

117 

378 

1.6 

1-Dec-ll 

17.2 

33.9 

1 

15-Dec-ll 

15.3 

38.8 

0.2 

15-Dec-ll 

16.4 

44.0 

0 

l-Jan-12 

29.2 

76.0 

8.8 

l-Jan-12 

166 

391 

2 

15-Jan-12 

31.2 

50.3 

11 

15-Jan-12 

31.5 

40.0 

11 

l-Feb-12 

42.4 

56.9 

24 

l-Feb-12 

39.9 

46.0 

21 

15-Feb-12 

30.5 

48.7 

10 

15-Feb-12 

29.7 

39.4 

9 

l-Mar-12 

40.4 

61.9 

20.0 

l-Mar-12 

34.5 

50.6 

14 

15-Mar-12 

33.3 

73.7 

13.0 

15-Mar-12 

33.4 

74.6 

13 

l-Apr-12 

32.2 

59.3 

12 

l-Apr-12 

38.2 

78.2 

16 

15-Apr-12 

10.0 

25.2 

0.00 

15-Apr-12 

9.2 

21.1 

0 

l-May-12 

25.7 

57.7 

5.2 

l-May-12 

20.8 

51.9 

2 

15-May-12 

14.7 

45.3 

0.1 

15-May-12 

12.2 

31.7 

0 

l-Jun-12 

18.1 

29.4 

1.1 

l-Jun-12 

22.3 

41.3 

3 

15-Jun-12 

11.0 

18.6 

0.00 

15-Jun-12 

10.7 

25.4 

0 

l-Jul-12 

7.2 

29.5 

0.00 

l-Jul-12 

7.6 

31.5 

0 

15-Jul-12 

10 

37 

2.6 

15-Jul-12 

12.5 

35.6 

0 
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Table  5.  YM  39  and  YM  38  Decibel  Values  Used  in  Figure  51 .  . 


YM  39 

YM  38 

Date 

dBatlHz 

dBat5  Hz 

Date 

dB  at  1  Hz 

dBat5  Hz 

13-Jun-08 

-135 

-154 

13-Jun-08 

-147 

-154 

14-Jun-08 

-158 

14-Jun-08 

-153 

l-Jul-09 

l-Jul-09 

-152 

15-Jul-08 

-86 

-103 

15-Jul-08 

-143 

-153 

l-Aug-08 

-97 

-113 

l-Aug-08 

-144 

-148 

15-Aug-08 

-96 

-110 

15-Aug-08 

-144 

-154 

l-Sep-08 

-93 

-108 

l-Sep-08 

-153 

15-Sep-08 

-95 

15-Sep-08 

-149 

l-Oct-08 

-93 

-111 

l-Oct-08 

-141 

-151 

15-Oct-08 

-100 

-116 

15-Oct-08 

-152 

l-Nov-08 

-96 

-111 

l-Nov-08 

-151 

15-Nov-08 

-99 

-113 

15-Nov-08 

-149 

13-Dec-08 

-122 

-137 

13-Dec-08 

-152 

14-Dec-08 

-122 

-126 

14-Dec-08 

-150 

l-Jan-09 

l-Jan-09 

-150 

15-Jan-09 

-114 

15-Jan-09 

-151 

l-Feb-09 

-128 

-151 

l-Feb-09 

-146 

-154 

15-Feb-09 

-143 

15-Feb-09 

-152 

l-Mar-09 

l-Mar-09 

-151 

15-Mar-09 

-122 

-141 

15-Mar-09 

-139 

-148 

l-Apr-09 

-123 

-143 

l-Apr-09 

-143 

-150 

15-Apr-09 

-130 

-132 

15-Apr-09 

-140 

-148 

l-May-09 

l-May-09 

-152 

15-May-09 

-148 

15-May-09 

-153 
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Table  6.  7D  J65A  and  7D  J73A  Used  in  Figure  52. 


7DJ65A 

7DJ73A 

Date 

dBatlHz 

dBat5  Hz 

Date 

dB  at  1  Hz 

dBat5  Hz 

1-Nov-ll 

-116 

-151 

1-Nov-ll 

-116 

-155 

15-Nov-ll 

-118 

-154 

-118 

-156 

1-Dec-ll 

-117 

-147 

-118 

-155 

15-Dec-ll 

-112 

-154 

15-Dec-ll 

-155 

l-Jan-12 

-114 

-152 

l-Jan-12 

-153 

15-Jan-12 

-115 

-152 

15-Jan-12 

-149 

l-Feb-12 

-115 

-148 

l-Feb-12 

-114 

-149 

15-Feb-12 

-118 

-150 

15-Feb-12 

-111 

-150 

l-Mar-12 

-114 

-153 

l-Mar-12 

-112 

-146 

15-Mar-12 

-118 

-150 

-109 

-149 

l-Apr-12 

-112 

-147 

-109 

-150 

15-Apr-12 

-121 

-148 

15-Apr-12 

-115 

-156 

l-May-12 

-114 

-153 

l-May-12 

-110 

-150 

15-May-12 

-118 

-150 

15-May-12 

-114 

-153 

l-Jun-12 

-124 

-153 

-117 

-148 

15-Jun-12 

-128 

-152 

-123 

-150 

l-Jul-12 

-121 

-149 

l-Jul-12 

-116 

-149 

15-Jul-12 

-119 

-148 

15-J  u  1-12 

-108 

-152 
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